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Report  Summary 


This  report  summarizes  activities  and  work  completed  between  September 
1,  1993  and  October  31,  1996.  In  addition  to  this  summary  it  contains  two 
parts  consisting  of  various  details. 

The  experimental  portion  of  this  project  is  supervised  by  O.  F.  Devereux 
and  conducted  by  P.  C.  Su.  As  indicated  in  the  earlier  reports,  it  is  de¬ 
voted  to  the  investigation  of  initiation  and  growth  kinetics  of  corrosion  in 
airframe  alloys,  including  joined  and  protectively  coated  specimens.  The 
goal  is  to  understand  fundamental  mechanisms  by  which  corrosion  damage 
develops  in  order  to  assist  in  design,  maintenance,  and  performance  evalua¬ 
tion  of  airframe  systems.  The  general  investigative  procedure  is  based  on  the 
principles  of  electrochemical  impedance  spectroscopy,  EIS.  These  principles 
together  with  the  problems  which  we  address  by  this  method  have  been  out¬ 
lined  in  the  earlier  reports  and  are  succinctly  summarized  in  the  introduction 
to  Part  1.  The  remainder  of  Part  1  is  devoted  to  a  detailed  description  of 
the  work  which  was  carried  out.  This  work  is  also  described  in  articles  7,  8, 
and  9  on  the  list  of  technical  reports. 

The  theoretical  portion  concerns  (i)  the  study  of  a  mathematical  model 
of  corrosion  conducted  by  G.  Hernandez  and  (ii)  investigation  of  certain  al¬ 
gorithmic  aspects  of  tomography  and  signal  and  image  processing  supervised 
by  W.  R.  Madych. 

The  objective  of  (i)  was  to  determine  the  theoretical  consequences  of 
a  well  known  mathematical  model  of  corrosion  in  alloys  and  the  details  of 
obtained  results  were  summarized  in  article  1  on  the  hst  which  was  included 
in  the  second  annual  report.  . 

The  work  associated  with  (ii)  is  detailed  in  articles  2,  3,  4,  and  5  on  the 
list.  Articles  2,  3,  and  4  were  summarized  and  included  in  the  first  two  annual 
reports.  Article  5  concerns  computed  tomography,  is  contained  in  Part  2  of 
this  report,  and  may  be  briefly  described  as  follows:  In  order  to  approximate 
a  small  portion  of  a  cross  section  conventional  reconstruction  algorithms 
require  the  averages  over  all  lines  which  pass  through  the  entire  cross  section. 
In  this  work  we  show  how  certain  high  frequency  components  of  a  portion 
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of  a  cross  section  can  be  approximated  using  only  averages  over  lines  which 
pass  through  or  close  to  the  region  of  interest.  The  resulting  high  frequency 
information  is  often  sufficient  for  the  detection  of  flaws,  abnormalities,  etc. 
Our  method  relies  on  the  notions  associated  with  the  so-called  continuous 
wavelet  transform  and  the  observation  that  Radon’s  classical  formula  is  a 
summability  formula  with  an  integrable  convolution  type  kernel.  What  is 
also  signiflcant  is  that  our  method  can  be  numerically  realized  via  simple 
modifications  of  well-established  convolution-backprojection  algorithms. 


Technical  Reports 


1.  G.  E.  Hernandez,  Internal-external  corrosion  in  alloys,  tech.  rep. 

2.  I.  Koltracht,  W.  R.  Madych,  and  K.  Marinelli,  Image  magnification 
using  families  of  fast  orthogonal  transforms,  Int  J.  Imaging  Systems 
and  Technology,  Vol.  7,  (1996),  170-179. 

3.  R.  A.  Lorentz  and  W.  R.  Madych,  Translation  and  dilation  invariant 
subspaces  and  multiresolution  analyses,  to  appear  in  Applied  and  Com¬ 
putational  Harmonic  Analysis. 

4.  W.  R.  Madych,  Finite  orthogonal  transforms  and  multiresolution  anal¬ 
yses  on  intervals,  to  appear  in  J.  Fourier  Analysis  and  Application. 

5.  W.  R.  Madych,  Reconstruction  and  wavelet  transforms  for  quasilocal 
tomography,  tech.  rep.  submitted  for  publication. 

6.  W.  R.  Madych,  Wavelets  and  their  applications,  a  guest  editorial  in 
Int.  J.  Imaging  Systems  and  Technology,  Vol.  7,  (1996),  149-151. 

7.  P.  C.  Su,  O.  F.  Devereux  and  W.  Madych,  Impedance  Imaging  For 
Prediction  and  Detection  of  Airframe  Corrosion,  in  Structural  Integrity 
in  Aging  Aircraft,  C.  I.  Chang  and  C.  T.  Sun,  eds.,  ASME  AD- Vol.  47, 
pp.  247-251, 
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8.  P.  C.  Su  and  0.  F.  Devereux,  AC  and  DC  Electrochemical  Studies  of  A 
Lanolin-Based  Inhibitor  on  2024  Aluminum,  submitted  for  publication 
in  Corrosion. 

9.  P.  C.  Su  and  O.  F.  Devereux,  Impedance  Studies  of  the  Corrosion  of 
Painted  and  Jointed  2024  Aluminum,  in  preparation. 

Articles  1,  2,  3,  4,  and  7  together  with  suitable  summaries  have  already 
been  submitted  with  earlier  technical  reports  to  the  AFOSR.  The  material 
in  articles  8  and  9  is  described  in  detail  in  the  enclosed  report  by  Devereux 
and  Su  in  Part  1.  Report  5  is  enclosed  in  its  entirety  in  Part  2. 


Presentations 


The  following  formal  presentations  where  made  by  the  investigators  as¬ 
sociated  with  this  grant: 

1.  W.  R.  Madych,  Impedance  Imaging  for  Prediction  and  Detection  of 
Airframe  Corrosion,  Second  USAF  Aging  Aircraft  Conference,  Okla¬ 
homa  City,  Oklahoma,  May,  1994.  Report  of  work  of  P.  C.  Su,  O.  F. 
Devereux  and  W.  R.  Madych.  Recorded  in  the  proceedings  edited  by 
C.  I.  Chang. 

2.  W.  R.  Madych,  Finite  orthogonal  transforms  and  multiresolution  anal¬ 
yses  on  intervals.  Tomography  Meeting  at  the  Mathematical  Research 
Institute  in  Oberwolfach,  Germany,  Sept.  1994.  Presentation  of  some 
of  the  results  in  article  4  listed  above. 

3.  W.  R.  Madych,  Image  magnification  using  families  of  fast  orthogonal 
transforms.  Wavelet  Meeting  at  the  Mathematical  Research  Institute 
in  Oberwolfach,  Germany,  July.  1995.  Presentation  of  some  of  the 
results  in  article  2  listed  above. 

4.  P.  C.  Su,  Impedance  Imaging  for  Prediction  and  Detection  of  Airframe 
Corrosion,  San  Francisco,  Nov.  12-17,  1995.  Presentation  of  results  in 
the  article  7  listed  above. 
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5.  O.  F.  Devereux,  Detailed  Modeling  of  AC  and  DC  Laboratory  Measure¬ 
ments  for  Corrosion  Prediction  at  Ambient  and  Elevated  Temperatures, 
NACE  Mid-America  Corrosion  Conference,  Chicago,  Oct.  1-4,  1995. 
Report  of  work  of  K.  S.  Yeum  and  O.  F.  Devereux. 

6.  O.  F.  Devereux,  Experimental  and  Theoretical  Aspects  of  Corrosion 
Detection  and  Prevention,  AFOSR  NDE  program  review,  Warner-Robins 
ALC,  Robins  AFB,  April  2,  1996.  Report  of  work  of  P.  C.  Su,  0.  F. 
Devereux  and  W.  R.  Madych. 

7.  W.  R.  Madych,  Pointwise  inversion  of  Radon  transforms.  Special  ses¬ 
sion  on  Tomography,  AMS  Meeting  in  Lawrenceville,  N.  J.,  Oct.  1996. 
Discussed  certain  consequences  of  some  of  the  results  in  article  5  listed 
above. 

We  mention  that  the  principal  investigator,  W.  R.  Madych,  was  a  guest 
editor  for  the  special  issue.  (Fall  1996)  of  the  Intermtional  Journal  of  Imaging 
Systems  and  Technology  devoted  to  wavelets  and  their  applications  in  image 
processing. 
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Part  1 
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INTRODUCTION 


The  Aging  Aircraft  Problem. 

The  aging  aircraft  problem,  that  of  maintaining  in  service  our  military  and  civilian 
fleets  for  periods  longer  than  their  projected  service  life,  has  been  recognized  for  some 
time.  A  broad  view  of  the  technological  issues  associated  with  this  problem  has  been 
given  recently  by  Chang,  including  understanding  the  corrosion  processes  that  occur  in 
aircraft,  the  detection  of  corrosion  damage  and  the  prevention  of  corrosion,  and  the 
understanding,  detection,  and  prevention  of  damage  due  to  mechanical  processes  such  as 
fatigue.'  These  may  be  intertvraied,  as  in  the  fatigue  failure  of  structures  weakened  by 
corrosion. 

The  mechanisms  by  which  corrosion  occurs,  or  “the  forms  of  corrosion,”  have  been 
categorized  in  various  ways;^  those  of  particular  interest  to  the  aircraft  problem  include 
general  corrosion,  galvanic  corrosion  and  crevice  corrosion.  All  corrosion  is 
electrochemical  in  nature  and  can  be  imderstood  in  terms  of  a  galvanic  cell.  In  general 
corrosion  the  structure  is  uniformly  attacked  by  an  oxidizing  environment  such  as  moist 
air.  In  galvanic  corrosion  dissimilar  metals  are  involved;  these  exhibit  different  affinities 
for  reducing  agents  in  their  environment  (e.g.,  in  moist  neutral  environments,  dissolved 
oxygen)  forcing  on  a  metal  to  be  cathodic  (site  of  reduction)  and  the  other  to  be  anodic 
(site  of  oxidation).  Unless  the  anode  forms  a  highly  protective  oxide  film  by  virtue  of  the 
corrosion  process,  it  is  attacked  and  dissolves  in  the  environment.  In  crevice  corrosion  a 
restricted  geometry,  such  as  a  riveted  joint,  creates  a  cell  by  virtue  of  the  different 
environments  within  and  without  the  crevice.  Within  the  crevice  the  dissolving  metal 
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hydrolyses  to  form  a  solid  oxide,  e.g.  Al(OH)3,  which  isolates  the  crevice,  permitting  the 
environment  within  the  interior  of  the  crevice  to  be  different  from  the  bulk  medium 
exterior  to  the  crevice.  The  act  of  hydrolyzing,  e.g.,  Al"^^  +  3H2O  =>  Al(OH)3  +  also 
acidifies  the  crevice,  enhancing  the  solubility  of  the  corrosion  product  and  preventing 
formation  of  a  protective  oxide  film,  while  the  net  positive  charge  in  this  region  attracts 
the  ubiquitous  chloride  ion.  Cl',  which  further  unstabilizes  the  protective  oxide  and 
accelerates  the  corrosion  process.  In  applying  these  mechanisms  to  aircraft  corrosion  we 
must  understand  that  aluminum  alloys  make  up  much  of  an  aircraft’s  structure,  that  a 
mature  aircraft  may  have,  literally,  thousands  of  pounds  of  brackish  water  wdthin  its 
structure,  that  aluminum  is  highly  reactive  and  is  stabilized  only  by  a  highly  protective 
AI2O3  film  which  in  turn  is  readily  attacked  by  Cl',  that  aircraft  structures  are  replete  with 
crevices,  such  as  those  formed  by  riveted  joints,  and  the  use  of  steel  fasteners  in 
juxtaposition  to  aluminum,  promoting  galvanic  corrosion,  is  commonplace.  The  detection 
of  skin  corrosion,  i.e.,  that  of  areas  of  sheet  aluminum,  is  being  addressed  by  ultrasonic, 
x-ray  diffraction  and  eddy  current  technology;  the  complex  geometry  of  joints  renders 
application  of  these  techniques  difficult  at  best,  useless  at  worst. 

Corrosion  of  Aluminum  Alloys. 

Aluminum  and  its  alloys  are  generally  somewhat  resistant  to  general  corrosion  in 
near-neutral  waters  due  to  the  highly  stoichiometric  and  thus  protective  oxide  film  that 
will  form  spontaneously  or  via  a  deliberate  process  such  as  anodizing  or  chemical 
conversion.  This  film  is,  however,  readily  penetrated  by  chloride  ion,  causing  pitting 
corrosion  which  proceeds  by  a  mechanism  somewhat  akin  to  crevice  corrosion.  Pitting  is 
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a  cosmetic  issue  in  aircraft,  but  may  also  become  a  structural  issue  through  the  creation  of 
stress  raisers.  Aluminum  is  susceptible  to  crevice  corrosion  because  of  the  acidification 
which  occurs  within  the  crevice  geometry  which  prevents  formation  of  protective  oxide, 
and  of  the  attraction  of  the  ubiquitous  and  detrimental  environmental  chloride  by  the 
metal  ions  within  the  crevice.  In  galvanic  couples  involving  brackish  waters  aluminum  is 
more  active,  or  anodic,  than  other  common  metals  except  for  magnesium  and  zinc,  but  its 
galvanic  behavior  is  not  always  predictable  due  to  the  extremely  passive  nature  of  its 
oxide  film.  Corrosion  of  aluminum  beneath  protective  organic  coatings  is  likely  to  occur 
in  a  manner  somewhat  analogous  to  crevice  corrosion.  Common  aircraft  alloys  such  as 
2024  and  7075  are  also  subject  to  exfoliation  corrosion,  detectable  by  ultrasonic  means, 
and  to  stress  corrosion  cracking.  The  latter  mode  of  failure  occurs  in  high  strength 
materials  in  particular  environments  when  the  material  is  stressed  above  a  critical  level;  it 
is  commonly  controlled  in  aircraft  by  assuring  that  service  loads  are  maintained  below  the 
critical  stress  level. 

Corrosion  Preventive  Inhibitor. 

Corrosion  preventive  inhibitors  are  frequently  used  for  protection 
of  aircraft  structures.  Current  corrosion  preventive  measures  rely  heavily  upon  inhibitor 
systems  which  use  volatile  solvent  carriers,  e.g.,  MIL-C-81309  in  which  the  active 
inhibitor  is  a  petroleum  sulfonate  and  which  forms  a  soft  “oil-like”  film  on  the  surface, 
and  MIL-C-85054  which  has  a  similar  base  but  also  contains  silicone  and  alkyd  resins 
and  cures  to  form  a  relatively  hard  film.  In  addition  to  the  environmental  liability  of  the 
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volatile  organic  carrier,  these  systems  are  not  suitable  for  moving  parts  and  afford  only 
temporary  protection,  yet  are  difficult  to  remove. 

Interest  in  new  protection  systems  for  aluminum  alloys  is  motivated  by  concern  for 
both  improved  protection  of  the  aircraft  and  for  the  environment.  In  this  work,  the  ability 
of  a  wool  wax-based  commercial  inhibitor.  Fluid  Film®  has  been  evaluated.  Wool  wax 
consists  primarily  of  lanolin,  whose  major  components  are  the  esters  RCOOR'  and 
R"COHCOOR'.  It  is  environmentally  benign,  and  is  applied  as  a  water-based  emulsion, 
thereby  avoiding  liability  of  volatile  organic  solvents. 

Electrochemical  Impedance  Spectroscopy  in  Corrosion  Research. 

The  corrosion  of  a  metal  is  an  electrochemical  process  and  involves  the  transport  of 
both  ions  and  electrons  in  the  vicinity  of  the  interface  between  the  metal  and  an 
electrolytic  environment.  It  has  a  resistive  aspect,  described  by  the  field  dependence  of 
the  anodic  and  cathodic  reactions,  and  a  capacitive  aspect  arising  from  a  non- 
homogeneous  distribution  of  ions  and  conduction  electrons  in  the  vicinity  of  the  interface. 
Therefore,  the  behavior  of  a  corroding  metal  can  often  be  represented  by  the  electrical 
behavior  of  a  simple  circuit  comprised  of  resistors,  capacitors  and  other  simple  elements 
such  as  shown  in  Figure  2  (a). 

In  Electrochemical  Impedance  Spectroscopy  (EIS)  a  small  AC  signal  (5-20  mv)  is 
applied  between  the  specimen,  termed  the  working  electrode,  and  an  inert 
counterelectrode  and  the  magnitude  and  phase  of  the  resulting  current  is  measured.  From 
this  information,  the  complex  impedance  of  the  interface  can  be  determined  over  a  range 
of  frequencies,  and  values  for  components  of  a  corresponding  equivalent  circuit  can  be 
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determined  through  fitting  the  impedance  spectrum  to  the  set  of  equations  describing  the 
model  circuit  via  the  use  of  software  such  as  that  developed  by  Boukamp/  Although  EIS 
may  be  used  to  quantify  corrosion  rates  over  a  broad  range  from  the  very  fast  to  the 
vanishingly  small,  it  is  especially  suited  to  systems  which  experience  low  corrosion  rates, 
i.e.,  well-protected  systems.  In  well-protected  systems  such  as  painted  or  anodized 
metals,  the  metal/electrolyte  interface  behaves  as  a  resistor,  Rp,  and  capacitor,  Cg,  in 
parallel,  and  the  electrolyte  behaves  as  a  series  resistor,  R*,  such  as  in  the  model  circuit 
shown  in  Figure  2  (b).  The  resistance  Rp  includes  the  corrosion  rate  of  the  metal  and 
electronic  leakage  through  the  protective  film,  while  the  capacitance  is  that  of  the 
protective  film  and  is  frequently  seen  to  reflect  the  diffusion  of  water  into  the  film.  As 
noted,  interfaces  are  characterized  by  a  non-uniform  distribution  of  mobile  charges,  or 
“electrical  double  layer,”  which  may  be  characterized  as  a  capacitor,  C<,i,  in  series  with 
the  film  but  may  not  be  seen  due  to  its  typically  large  magnitude  (capacitors  in  series  add 
as  reciprocals). 

Degradation  of  protective  coatings  and  corrosion  of  metal  substrates  can  be  assessed 
from  changes  of  the  values  of  equivalent  circuit  parameters  determined  through  EIS.  One 
advantage  of  EIS  over  traditional  DC  polarization  is  its  non-destructive  nature;  the  use  of 
a  small  excitation  signal  permits  determination  of  long-term  corrosion  properties  without 
perturbing  the  system.  Current  computer-driven  instrumentation  renders  this  technique 
both  user-friendly  and  rapid,  and  permits  explicit  evaluation  of  the  individual  electrical 
characteristics  of  even  highly  resistive  interfaces  such  as  well-protected  airframe 
structures. 
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Research  Plan. 


This  program  has  focused  upon  crevice  corrosion  in  simulated  aircraft  structures  as  a 
critical  mode  of  failure  that  is  not  adequately  addressed  by  the  more  traditional  NDT 
methods  now  under  investigation.  In  addition,  it  has  employed  electrochemical  methods 
in  a  study  of  the  efficacy  of  lanolin-based  inhibitors,  with  particular  attention  to  crevice 
and  to  previously  corroded  structures. 
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EXPERIMENTAL  OUTLINE 


The  aluminum  alloy  2024-T3  with  and  without  a  MIL-P-23377+MIL-C-85285  epoxy 
primer  and  polyurethane  topcoat  system  were  being  used  in  this  study.  These  coated 
specimens  were  provided  by  the  Navy  Air  Warfare  Center,  Patuxent  River,  Maryland, 
and  were  designed  to  meet  the  extremely  corrosive  service  environment  of  Naval  aircraft. 
Bare  2024-T3  aluminum  specimens  were  polished  to  a  600  grit  finish,  washed  with 
double  distilled  water  and  degreased  with  acetone.  Inhibitor-coated  specimens  were 
prepared  by  applying  the  wool  wax-based  inhibitor.  Fluid  Film®,  by  brush,  providing  a 
weight  gain  after  drying  of  ~9.5  mg/cm^. 

The  impedance  cell  was  a  2  in.  dia.,  4  in.  high  (5  cm  x  10  cm)  Pyrex®  cylinder, 
sealed  to  the  3"x3"(7.5  cm  x  7.5  cm)  x  1.5  mm  specimen  via  a  Viton®  0-ring,  which 
provided  an  exposed  specimen  area  of  18.3  cm  .  The  cell  contained  a  simulated  service 
environment  comprised  of  150  cc  of  0.5M  NaCl  in  double  distilled  water.  A  crevice  was 
simulated  by  affixing  a  Teflon®  cap  screw  with  a  shoulder  through  the  exposed  surface 
of  the  specimen  using  a  predetermined  torque  (160  inch  ounces).  The  inhibitor-coated 
specimens  were  air-dried  for  three  hours  before  assembling  into  the  cell,  which  was  open 
to  the  atmosphere  to  permit  natural  aeration.  The  aluminum  specimen  comprised  the 
working  electrode,  the  coimterelectrode  was  a  316L  stainless  steel  plate  and  the  reference 
electrode  was  a  saturated  calomel  electrode.  Impedance  spectra  were  obtained  with  a 
Solatron  Model  1250  Frequency  Response  Analyzer  coupled  with  a  Solatron  Model  1286 
Potenstiostat  and  an  on-line  computer.  Z-plot®  software  was  used  to  control  the 
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experiments,  and  equivalent  circuit  modeling  was  performed  with  Boukamp’s 
EQUICRT  software.^ 

Measurement  of  the  open  circuit  potential  and  DC  electrode  polarization  studies  were 
also  conducted  in  a  glass  cell  with  an  air-purged  0.5M  NaCl  aqueous  solution  as  the 
electrolyte,  a  1  cm  2024-T3  disk  as  working  electrode,  a  saturated  calomel  electrode  as 
reference  electrode,  and  a  graphite  rod  as  the  counterelectrode.  These  experiments  were 
performed  with  a  Princeton  Applied  Research  Potentiostat/Galvanostat  Model  273 
interfaced  to  a  personal  computer,  using  PAR  352  Corrosion  Analysis  Software  and  a 
polarization  scan  rate  of  0.17mv/s. 
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PROGRAM  ACCOMPLISHMENTS 


A.  Outline  of  Achievements. 

The  work  performed  imder  this  contract  consisted  of  a  number  of  tasks: 

1.  Define  the  fundamental  impedance  response  of  the  aluminum  airframe  alloys  by 
performing  baseline  EIS  measurements  in  saline  environment;  and  determine  the  time 
dependence  of  corrosion  in  systems  of  interest  via  dynamic  impedance  analysis.  This 
task  has  been  completed. 

2.  Apply  these  methods  to  specific  structures  representative  of  airframe  geometries,  i.e., 
joined  and  protected  coated  specimens.  Aluminum  specimens  with  simulated  joint 
were  studied  in  corrosive  environment. 

3.  Develop  a  methodology  to  distinguish  between  a  corroding  structure  and  an  inert 
stmcture,  e.g.,  between  a  sound  riveted  joint  and  a  riveted  joint  undergoing  hidden 
corrosive  attack.  The  extent  to  which  the  different  structures  have  corroded  can 
be  determined  from  the  corresponding  type  of  equivalent  circuit,  and  from  the 
magnitude  of  the  equivalent  circuit  parameters. 

4.  Evaluate  the  effectiveness  of  a  new  eviromentally  innocuous  commercial  inhibitor, 
“Fluid  Film”,  with  electrochemical  techniques  as  an  alternative  for  the  protection  of 
the  aircraft  structures  from  corrosion.  DC  and  AC  electrochemical  techniques  have 
been  applied;  results  indicate  that  the  inhibitor  can  effectively  protect  aluminum 
from  corrosion  even  if  initial  corrosion  has  already  occurred. 

These  results  are  described  in  the  following  pages: 
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B. 


Detail  Tasks  Summaries: 


1.  Baseline  BIS  response  and  time  dependence  of  EIS  behavior  for  2024-T3  airframe 

aluminum  allov  in  saline  environment. 

Extensive  EIS  data  have  been  collected  on  the  aluminum  alloy  2024-T3  immersed 
in  aqueous  0.5M  NaCl  in  the  (a)  as  received,  (b)  painted,  and  (c)  treated  with  a 
commercial  wool  waxed-based  inhibitor  conditions.  These  impedance  spectra  are  well 
behaved  and  reproducible,  and  can  be  fit  with,  a  high  degree  of  accuracy  to  appropriate 
equivalent  circuits.  Typical  Nyquist  and  Bode  representations  for  bare  aluminum  are 
shown  in  Figure  1  (a)  and  (b);  they  are  modeled  with  the  Randles  equivalent  circuit, 
shown  in  Figure  2  (a),  and  the  dependence  of  the  equivalent  circuit  parameters  with 
immersion  is  given  in  Table  1.  In  Randles  circuit  R^  represents  a  series  or  electrolyte 
resistance,  Cji  is  the  capacitance  of  the  double  layer,  R^  is  the  charge  transfer  resistance, 
and  W  is  the  Warburg  impedance,  given  by  ^ 

W=(<y/a)‘'^)(l-j)tanh[;i(j®/D)*^] 


where  co  is  the  circular  frequency,  a  is  the  “Warburg  coefficient,”  h  is  the  diffusion 
thickness  and  D  is  the  relevant  diffusivity. 

Figure  3  (a)  and  (b)  display  impedance  spectra  for  painted  aluminum  with  an  intact 
coating  and  damaged  coating,  respectively,  and  model  curves  using  with  the  equivalent 
circuit  shown  in  Figure  2  (b).  These  spectra  display  a  effect  of  the  presence  of  flaws  in 
the  coating  as  shown  in  Figure  3  (b).  Such  behavior  may  or  may  not  be  accompanied  by  a 
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visible  blister  or  other  artifact  in  the  paint  and  in  long  term  exposure  tests  is  generally 
discernible  via  EIS  prior  to  the  appearance  of  any  visible  defect.  The  differences  in  these 
intact  and  damaged  coatings  is  shown  quantitatively  via  the  effect  on  model  circuit 
parameters  in  Figures  4  through  7.  The  coating  capacitance,  C^,  increases  markedly 
during  the  first  30  days  of  exposure,  Figure  4;  this  behavior  is  consistent  with  water 
uptake  by  the  coating  and  is  similar  for  the  two  specimens  because  this  phenomenon  is 
independent  of  damage.  The  step-like  appearance  of  the  capacitance  is  attributed  to  “non- 
Fickian”  diffusion  of  water  into  the  film  caused  by  successive  occupation  of  different 
types  of  site  in  the  film  by  water.  Figure  5  compares  the  film,  or  “pore,”  resistance  of  the 
damaged  and  intact  paint,  a  parameter  which  reflects  the  DC  conductivity  of  the  paint. 
For  the  intact  film  this  decrease  in  accord  with  the  water  adsorption  shown  in  Figure  5; 
the  damaged  film  displays  a  significantly  lower  resistance  that  rapidly  deteriorates  as 
major  film  flaws  develop.  Both  parameters,  and  Rpo,  thus  quantitatively  indicate 
degradation  of  protective  coatings  by  physical  damage,  intrinsic  defects  or  diffusive 
invasion.  Figure  6  compares  the  polarization  resistance  of  the  metal  substrate  for  the  two 
coatings.  This  parameter  is  a  direct  measure  of  the  resistance  of  the  system  to  corrosion  in 
that  its  reciprocal,  1/Rp,  is  representative  of  corrosion  rate.  Polarization  resistance 
gradually  deteriorates  with  time  in  both  systems,  suggesting  continued  degradation  after 
the  coatings  become  saturated  with  water;  and  Rp  is  significantly  lower  for  the  damaged 
system,  indicative  of  very  little  protection.  The  variation  in  the  double  layer  capacitance 
with  exposure  time  shown  in  Figure  7  is  marked.  Initially  very  low,  the  value  of  Cji  for 
the  intact  coating  increases  by  nearly  three  orders  of  magnitude  to  approximately  5 
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|aF/cm^.  The  increase  in  this  parameter  with  time  indicates  the  development  of  an  ionic 
space  charge  within  the  paint  film;  the  much  larger  values  seen  for  the  damaged  film 
represents  actual  contact  of  water  with  the  metal. 
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2.  ETS  study  of  jointed  and  protected  aluminum  airframe  structures. 


Extensive  testing  similar  to  that  performed  on  simple  specimens  has  been  performed 
on  a  structural  model  comprised  of  the  painted  A1  2024-T3  vsith  a  simulated  rivet,  i.e.,  a 
stepped  nylon  cap  screw  affixed  to  the  test  plate  with  a  predetermined  torque,  to  provide  a 
reproducible,  defined  crevice  geometry  without  introducing  the  artifact  of  a  galvanic 
couple.  Figure  8  (a)  and  (b)  depict  Nyquist  and  Bode  representation  of  the  EIS  of  such  a 
joint  over  a  month  of  immersion  in  0.5M  NaCl.  The  short  term  behavior  is  similar  to  that 
of  the  simple  specimen  of  Figure  3  vdth  two  time  constants  displayed,  but  upon  exposure 
the  impedance  decreases  indicating  a  region  of  weaker  protection  associated  with  the 
crevice  which  appears,  however,  to  be  stable. 

With  a  riveted  joint  displaying  visible  corrosion,  shown  in  Figure  9  (a)  and  (b),  we 
see  a  marked  reduction  in  impedance  magnitude,  but  also  a  characteristic  third  time 
constant  which  appears  to  represent  a  region  of  active  corrosion  within  the  crevice.  The 
equivalent  circuit  of  Figure  2  (c)  was  used  to  model  the  corroding  crevice.  Figure  10 
illustrates  the  comparison  of  corroding  and  non-corroding  joints;  for  the  corroding 
crevice,  the  impedance  value  after  4  days  immersion  in  the  test  electrolyte  is  already 
much  smaller  than  that  of  the  non-corroding  crevice  after  36  days  immersion.  Circuit 
parameters  for  corresponding  to  the  corroding  joint  in  Figure  9  are  shown  in  Figure  1 1 
and  Figure  12.  The  coating  capacitance,  Cj.,  and  the  double  layer  capacitance,  Cji,  both 
increase  with  time,  this  can  be  attributed  to  the  water  diffusion  into  the  coating  as  with 
case  of  the  coated  sample.  The  crevice-related  capacitance,  C^y,  increases  with  time  up  to 
30  days  then  decreases.  The  capacitance  values  associated  with  the  corroding  crevice  are 
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abnormally  high  (~niF/cm^).  Such  large  capacitance  values  have  been  interpreted  as 
“supercapacitance”,’  and  may  be  associated  with  Faradaic  (or  corrosive)  processes  at  the 
surface  which  mimic  the  charging  characteristics  of  a  capacitor,  or  as  a 

o 

“pseudocapacitance”  which  may  be  associated  with  an  electrosorption  process.  The  pore 
resistance,  Rpo,  and  the  polarization  resistance,  Rp,  both  decrease  with  immersion  time  as 
with  the  simple  coated  specimens.  The  crevice-related  resistance,  R^v,  may  be  related  to 
the  charge  transfer  resistance  of  the  crevice  and  displays  relatively  small  values  which 
decrease  with  immersion  time,  suggesting  minimum  corrosion  protection.  It  has  been 
reported  that  the  pH  value  of  the  fluid  within  the  crevice  for  neutral  dilute  sodium 
chloride  solution  is  near  2,’  highly  solubilizing  to  the  protective  oxide  layer  on 
aluminum.  Tests  of  macro  specimens  in  similar  solutions  confirm  the  role  of  the  crevice 
environment  in  facilitating  corrosion.  The  equivalent  circuit  parameters  for  bare  2024-T3 
in  0.5M  NaCl  aqueous  solution  pH  2  are  shown  in  Table  2,  Comparing  R^j  values  in 
Table  1  with  those  in  Table  2,  the  charge  transfer  resistance  in  the  acidic  environment  is 
at  least  10  times  smaller  than  that  in  neutral  solution.  Thereby,  the  EIS  technique  allows 
us  to  discern  the  extent  of  corrosion  in  joined  and  protected  aluminum  structures,  and  the 
presence  of  an  actively  corroding  crevice  can  be  clearly  detected  in  contrast  to  other 
techniques. 
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3.  Methodology  to  differentiate  a  corroded  structure  from  a  sound  structure 


In  the  previous  task  1,  we  showed  that  the  degree  of  protection  afforded  by  a  coating 
on  an  aliuninimi  structure  can  be  quantified  by  the  EIS  technique  in  terms  of  the 
magnitude  of  equivalent  circuit  parameters  such  as  Cc  and  Rc,  before  any  visual 
degradation  is  appeared.  Through  task  2,  hidden  corrosion  such  as  crevice  corrosion  in  a 
fabricated  structure  was  shown  to  be  detected  by  the  analysis  of  equivalent  circuit 
parameters,  Ccv  and  R^v  obtained  via  EIS.  Thereby,  if  equivalent  circuit  analysis  are 
obtained  with  EIS  data  from  an  in-service  structure,  information  about  the  condition  of 
the  structure  with  respect  to  existing  corrosive  damage  can  be  easily  determined. 
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4.  Apply  electrochemical  techniques  to  evaluate  the  effectiveness  of  an  environmentallY 
innocuous  corrosion  inhibitor  on  airframe  aluminum  alloy, 

The  variation  of  open  circuit  potential  (OCP)  with  immersion  time  for  as-polished 
and  inhibitor-coated  specimens  is  shown  in  Figure  13.  The  OCP  of  the  bare  specimen 
drops  from  the  initial  value  of -6 10  mV  SCE  and  approaches  a  nearly  steady  value  of 
-740  mV  SCE  after  ~20  h.  The  bare  sample  specimen  shows  considerable  fluctuation  in 
OCP  value  between  5  and  15  hours  of  immersion,  suggestive  of  vigorous  initial  pitting 
activity  that  latter  stabilizes.’  The  initial  decrease  of  the  OCP  of  the  bare  sample  is 
consistent  with  an  increase  of  the  average  anodic  current  due  to  rapid  localized 
dissolution  rates  at  pits  sites,  with  the  subsequent  steady  value  indicative  of  a  stable  pit 
structure.  The  OCP  of  the  inhibitor-coated  aluminum  increases  monotonically  with 
immersion  time  and  approaches  a  steady  state,  indicative  of  an  absence  of  pitting  and 
suggestive  of  a  diminishing  anodic  reaction  rate.  This  would  be  expected  if  the  inhibitor 
film  encourages  a  slow  but  uniform,  i.e.,  pit-free,  growth  of  aluminum  oxide. 

Polarization  curves  for  bare  2024-T3  aluminum  in  aerated  0.5M 
aqueous  NaCl  and  for  inhibitor-coated  2024-T3  aluminum  in  the  same  environment  are 
displayed  in  Figure  14.  The  corrosion  potential  of  bare  aluminum  is  higher  than  that  of 
the  inhibitor-coated  aluminum,  and  corresponds  to  the  OCP  upon  initial  inunersion  in 
Figure  13.  The  current  densities  for  bare  aluminum  are  about  4  orders  of  magnitude  larger 
than  those  for  coated  aluminum,  indicative  of  the  protective  nature  of  this  inhibitor. 

Polarization  data  for  the  bare  and  the  coated  aluminum  electrodes  were  analyzed 
via  the  modeling  technique  of  Yeum  and  Devereux,*’* "  and  are  shown  in  Figures  15  and 
16,  respectively,  together  with  model  curves.  In  this  model  the  polarization  potential  is 
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the  sum  of  the  potentials  attributable  to  electron  transfer  (activation)  polarization, 
diffusion  (concentration  polarization),  and  ohmic  drop  in  the  electrolyte.  For  the  jth 
reaction  this  is  expressed  as: 

ij  =  Sj(l-Sjij/i'j)'"jexp[Sj(V-V*j)/bj]exp(-iARj/bj) 

where  Sj  is  +1  for  an  anodic  reaction  and  -1  for  a  cathodic  reaction,  i^j  is  the  limiting 
diffusion  current,  mj  is  the  reaction  order,  2.303bj  is  the  Tafel  slope,  A  is  the  surface  area 

]|c 

of  the  electrode  and  R  is  the  cell  resistance.  The  defined  parameter  V  j  combines  the 
functions  of  the  equilibrium  potential  and  the  exchange  current  and  is  given  by 


V*=  V«f  -  Sjbjln(Zj3kj  )-Sjbjln(  H  ) 

i=1 

where  V^gf  represents  the  potential  of  the  reference  electrode  relative  to  the  electrolyte,  Zj 
is  the  munber  of  electrons  transferred  in  the  rate  controlling  step,  3  is  the  Faraday,  kj  is 
the  chemical  rate  constant,  and  the  aj  are  the  activities  of  the  reacting  species.  The  total 

current  experienced  by  the  electrode  is  the  sum  of  the  currents  carried  by  the  individual 
reactions: 


nj 

i  =  l  Eh 

j=l 


The  technique  utilizes  a  nonlinear  regression  procedure  to  optimize  the  values  of  the 
model  parameters. 

The  polarization  curve  for  bare  2024-T3  aluminum  in  aerated  0.5M  NaCl,  Figure 
15,  is  modeled  by  three  reactions:  C2,  reduction  of  hydrogen  ion,  modeled  by  electron 
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transfer  polari2ation  and  ohmic  drop;  Cl,  reduction  of  dissolved  oxygen,  modeled  by 
electron  transfer  polarization  and  diffusion;  and  Al,  anodic  oxidation  of  aluminum, 
modeled  by  electron  transfer  polarization  and  ohmic  drop.  Although  the  model  can 
incorporate  passivation  there  was  no  evidence  of  passive  behavior,  perhaps  due  to  the 
almost  immediate  onset  of  chloride-induced  pitting.  Model  parameters  for  these  reactions 
are  listed  in  Table  3;  the  corrosion  current,  from  the  model,  is  0.104  mA/cm^. 

The  polarization  curve  for  inhibitor-coated  2024-T3  aluminum  in  aerated  0.5M 
NaCl,  Figure  16,  is  modeled  by  two  reactions:  Cl,  reduction  of  dissolved  oxygen,  and 
Al,  anodic  oxidation  of  aluminum,  both  controlled  by  electron  transfer  polarization  and 
diffusion.  Model  parameters  for  these  reactions  are  listed  in  Table  4;  the  corrosion 
current,  from  the  model,  is  9.35  nA/cm^,  four  orders  of  magnitude  smaller  than  that  of 
bare  aluminum. 

Impedance  spectra  were  recorded  after  various  times  of  immersion  in  naturally 
aerated,  aqueous  0.5M  NaCl;  Nyquist  and  typical  Bode  plots  are  shown  for  inhibitor- 
coated  specimens  in  Figure  17  (a)  and  (b).  Historically,  various  equivalent  circuit  models 
have  been  used  to  model  impedance  spectra  of  inhibitor-coated  aluminum,*^’  e.g.,  two 
RC  elements  in  series  or  in  parallel,  with  and  without  Warburg  elements,  euid  three  RC 
elements  in  parallel.  For  the  spectra  in  Figures  17  the  authors  adopted  the  equivalent 
circuit  model  shown  in  Figure  2  (d).  Although  these  spectra  appear  to  represent  but  two 
time  constants  (i.e.,  there  are  two  semicircular  elements),  two  are  needed  to  accurately 
model  the  depressed  higher  frequency  semicircle,  with  the  third  representing  the  low 
frequency  arc.  In  this  circuit  R^,  C<i|,  and  R^t  retain  their  significance  in  the  Randles 
circuit;  C,.  represents  the  capacitance  of  the  inhibitor  film,  Rp^  is  the  “pore  resistance”  of 
the  film;  and  Qx  and  R<,x  are  the  capacitance  and  resistance,  respectively,  of  the  oxide 
layer.  To  improve  the  fit  of  model  to  data,  in  all  cases  capacitors  were  replaced  by 
constant  phase  elements  (CPE)  of  impedance  ^ 
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^CPE  ~  (j®) 

where  a  is  a  constant  (-1  <  a  <  1)  and  Aq'*  corresponds  to  the  capacitance  if  a  is  equal  to 
1.  The  dependence  of  the  circuit  parameters  for  the  coated  specimens  (excepting  R^, 
which  is  not  significant)  with  immersion  time  is  shown  in  Figures  1 8  and  19. 

Figure  20  (a)  shows  a  Nyquist  representation  of  impedance  spectra  for  a  specimen 
that  had  been  immersed  in  the  brine  environment  for  8  days  prior  to  treatment  with  the 
inhibitor;  Figtire  20  (b)  shows  a  typical  Bode  plot  for  this  specimen.  The  equivalent 
circuit  model  of  Figure  2  (d)  was  used  to  fit  these  spectra;  the  corresponding  values  of  the 
equivalent  circuit  parameters  are  exhibited  in  Figure  21  and  22.  The  values  for  the 
coating  capacitance  and  the  double  layer  capacitance,  C5.  and  Cji,  for  this  “pre-corroded” 
specimen  are  larger  than  those  for  the  as-received  aluminum  treated  with  inhibitor,  while 
the  capacitance  of  the  oxide  film,  C^x,  is  smaller  for  the  specimen  treated  after  initial 
corrosion.  The  resistance  values  for  this  specimen.  Figure  22,  as  with  the  other  impedance 
specimens  display  somewhat  more  erratic  behavior  than  the  capacitance  values. 

Equivalent  circuit  in  Figure  2  (b)  describes  the  impedance  behavior  of  the 
inhibitor-coated  aluminum  substrate.  This  model  is  composed  of  three  distinct  RC 

elements,  whose  correlation  with  physical  features  of  the  system  is  readily  accomplished 
by  examining  the  magnitudes  of  the  measured  capacitance  values.  The  capacitance  is 

of  the  order  of  10'*°  F/cm^  and  remains  essentially  constant  throughout  the  experimental 
period..  We  may  assume  a  dielectric  constant  of  10  for  the  rather  polar  species  (waxes, 
alcohols,  and  free  fat  acids)  comprising  wool  wax  and  use,  for  the  parallel  plate  capacitor 
constituting  the  inhibitor  film. 


C  =  eso/5 

where  e  is  the  dielectric  constant  of  the  film,  Eq  is  the  permittivity  of  free  space 
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(8.85x10’^^  farads/meter),  and  5  is  the  film  thickness.  Using  the  measured  capacitance 
value,  the  calculated  film  thickness  is  of  the  order  of  0.01  cm,  in  excellent 
correspondence  to  the  mass  of  the  applied  inhibitor  film,  9.5  mg  cm"^.  This  film  is 
applied  as  an  aqueous  emulsion  and  air  dried  and,  as  the  value  of  C^.  remains  constant 
during  immersion,  the  film  does  not  appear  to  absorb  additional  water  during  immersion. 
The  value  of  C^x,  associated  with  the  capacitance  of  the  oxide,  ranges  fi-om  1.5  to  2  x  10’^ 

F/cm*^;  this  corresponds  to  an  oxide  film  thickness  of  4  to  5  nm  if  a  dielectric  constant  of 
9  is  assumed.  This  is  just  slightly  thicker  than  the  thickness  of  the  native  AI2O3  film,  2  - 

4  nm,  reported  by  Scully  et  al,'"'  suggesting  that  this  circuit  element  is  indeed  the 
capacitance  of  the  oxide  film  and  that  the  rate  of  growth  of  the  oxide  is  very  slow  for 
inhibitor-coated  aluminum.  Indeed,  in  contrast  to  the  vmtreated  alloy  the  inhibitor-coated 
specimens  did  not  exhibit  visible  corrosion  during  the  period  of  immersion.  We  assoeiate 
the  capacitance  C^i  with  the  double  layer,  although  the  value  of  Cqx  more  closely 
resembles  typical  double  layer  values  in  aqueous  systems.  However,  the  correlation  of  Cf 
with  the  wool  wax  layer  is  definitive,  and  the  value  of  Cjn  cannot  represent  the  oxide  film 
because  the  two  order  of  magnitude  reduction  in  thickness  implied  by  the  increase  in  this 
parameter  is  not  plausible.  Rather,  we  suggest  that  the  value  of  Cji,  inereasing  from  2x10"* 
F/cm  to  0.02  F/cm  during  the  month  of  immersion,  represents  the  diffusion  of 
electrolyte  into  the  wool  wax  film  to  create  a  very  diffuse  Gouy  layer  of  charge. 

As  noted,  the  values  of  the  resistive  elements  exhibit  a  somewhat  more  erratic 
behavior  than  do  the  capacitance  values  for  these  specimens.  The  oxide  resistance,  R^x, 
shows  an  overall  increase  which  would  be  associated  with  growth  of  aluminum  oxide 
were  it  not  for  the  invariance  of  the  corresponding  capacitance  value.  An  alternative 
explanation  is  an  increase  in  the  average  resistivity  of  the  oxide  film  through  selected 
closure  of  film  defects.  Indeed,  such  a  mechanism  would  be  an  explanation  for  the 
efficacy  of  this  inhibitor.  The  resistance  of  the  inhibitor  film  shows  an  initial  increase  but 
is  subsequently  constant,  suggesting  a  similar  defect  closure  process.  Aside  from  an 
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initial  order  of  magnitude  increase  in  value,  the  charge  transfer  resistance  mimics  the 
shape  of  the  oxide  resistance.  We  believe  it  plausible  that  this  is  not  a  real  effect  but  an 
artifact  of  the  modeling  process. 
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roNruisioN  of  major  findings 

•  The  EIS  behavior  of  aluminum  specimens  -  bare,  painted,  inhibited,  and  with  crevice 
geometry  -  in  brackish  water  is  well  behaved  and  suitable  for  modeling. 

•  Behavior  such  as  adsorption  of  water  by  oxide,  paint  or  inhibitor  films  is  clearly 
discernible  via  the  EIS  technique. 

•  EIS  model  parameters  are  markedly  different  between  damaged  and  undamaged 
protective  coatings. 

•  EIS  and  DC  electrode  polarization  studies  confirm  the  empirically  observed  benefits 
of  the  lanolin-based  inhibitors  on  airframe  aluminum  alloys. 

•  The  presence  of  a  corroding  crevice  on  the  test  specimen  is  readily  apparent,  both 
qualitatively  (three  time  constants)  and  quantitatively  (“supercapacitance”  behavior, 
low  polarization  resistance). 
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Table  1 

Equivalent  circuit  parameters  for  bare  2024-T3  aluminum. 


Immersion  time 
(day) 

Cdl 

(nF/cm^) 

n 

Ret 

(ohm  cm2) 

(ohm) 

h/D'/2 

(sec'A) 

■191 

1 

38.1 

0.923 

9140 

310.6 

15.14 

4700 

3 

39.4 

0.954 

11200 

346.0 

13.24 

4580 

5 

45.3 

0.948 

11600 

454.5 

13.55 

6160 

Table  2 

Equivalent  circuit  parameters  for  bare  2024-T3  aluminum  in  0.5M  NaCl  aqueouss 

solution  with  pH  =  2. 


n 

Ret  (ohm  cm"^) 

154.4 

0.885 

833 
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Table  3 

Polarization  curve  parameters  for  bare  2024-T3  aluminum. 


Reaction 

s 

V*  (mV) 

b(mV) 

m 

i®  (mA) 

R(ohm) 

1(C2) 

-1 

-858.6 

39.1 

- 

- 

1.18 

2  (Cl) 

-1 

-789.0 

97.1 

1 

0.236 

- 

3(A1) 

1 

-540.3 

12.7 

- 

11.29 

Table  4 

Polarization  curve  parameters  for  inhibitor-coated  2024-T3  aluminum. 


Reaction 

s 

V*  (mV) 

b(mV) 

m 

i®  (^A) 

R(ohm) 

1(C1) 

-1 

-1095.47 

130.93 

1 

0.0476 

“ 

2(A1) 

1 

-648.4 

100.95 

1 

- 
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Impedonce  (ohm*cm''2)  -Zimag  (ohm*cm''2) 


Figure  1  (a).  Experimental  and  model  Nyquist  plots  for  bare  2024-T3 
aluminum  in  aqueous  0.5M  NaCI  solution. 


Frequency  (Hz) 

Figure  1  (b).  Typical  Bode  plots  corresponding  to  Figure  I  (a). 


Phase  Angle  (degree) 


Figure  2  (a).  Equivalent  electrical  circuit  for  metal 
without  coating. 


Figure  2  (b).  Equivalent  electrical  circuit  for  metal  with 
coating. 


Figure  2  (c).  Equivalent  electrical  circuit  for  coated  metal 
with  simulated  joint. 


Figure  2  (d).  Equivalent  electrical  circuit  for  inhibitor- 
coated  metal. 
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Figure  14.  Experimental  curves  for  bare  and  inhibitor-coated  2024- r3  Figure  15.  Experimental  and  model  polarization  curves,  showing 

aluminum  in  aqueous  0.5M  NaCI.  individual  contributions  of  electrode  reactions,  for  bare  2024-13  aluminum 

in  aqueous  0.5M  NaCI. 
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Figure  18.  Dependence  of  capacitance  elements  for  Figure  17  (a) 

Figure  17  (b).  Typical  Bode  plot  corresponding  to  Figure  1 7  (a).  specimen  on  time  of  immersion  in  aqueous  0.5M  NaCI. 


R  (ohm*cm^2)  C  (F/cm^2) 


Figure  21.  Dependence  of  capacitance  elements  for  Figure  20  (a) 
specimen  on  time  of  immersion  in  aqueous  0.5M  NaCI  after  inhibitor 
application. 


Figure  22.  Dependence  of  resistance  elements  for  Figure  20  (a)  specimen 
on  time  of  immersion  in  aqueous  0.5M  NaCI  after  inhibitor  application. 
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Reconstruction  and  Wavelet  Transforms  for 
Quasilocal  Tomography 

W.  R.  Madych  * 


Abstract 

The  terms  “local”  and,  more  recently,  “pseudolocal  tomography” 
have  been  introduced  by  several  authors  and  refer  to  paradigms  for 
recovering  high  frequency  information  concerning  planar  densities  /  at 
points  X  from  data  which  consists  of  averages  of  /  over  lines  through 
X  or  close  to  x.  The  wavelet  transform  of  /  is,  in  some  sense,  an  ideal 
tool  for  analyzing  local  frequency  content.  In  this  paper  we  indicate, 
among  other  observations,  how  reconstruction  and  wavelet  transforms 
of  /  can  be  computed  from  such  data.  Furthermore  we  show  how 
such  transforms  can  be  numerically  realized  via  simple  modifications 
of  well-established  convolution-backprojection  algorithms.  As  part  of 
our  development  we  review  the  concepts  of  “local”  and  “pseudolocal 
tomography”  and  show  that,  in  effect,  these  notions  basically  involve 
the  computation  of  a  wavelet  transform.  The  results  in  this  paper  are 
based  on  the  observation  that  Radon’s  classical  inversion  formula  is  a 
summability  formula  with  an  integrable  convolution  type  summability 
kernel. 


1  Introduction 

1.1  Background 

The  values  of  a  sufficiently  well  behaved  scalar  valued  density  or  function 
/  defined  on  the  plane,  x  E  IR^,  can  be  recovered  from  values  of  its  Radon 
transform  via  Radon’s  celebrated  classical  formula 


*Department  of  Mathematics;  University  of  Connecticut;  Storrs,  CT  06269-3009.  Par¬ 
tially  supported  by  AFOSR  Grant  F49620-93-1-0565 
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where  Fx{q)  is  the  average  of  all  the  integrals  of  /  over  lines  which  are  distance 
q  from  x.  For  example^  see  [29,  page  245]  also  [4,  Appendix  A],  [8],  [22],  and 
[25]. 

It  is  well  known  that  formula  (1)  contains  features  which  make  it  difficult 
to  use  in  practical  applications.  Among  these  are  (i)  the  fact  that  the  integral 
in  the  formula  is  not  proper  and  (ii)  the  fact  that  the  inversion  is  not  local, 
namely,  the  evaluation  of  f{x)  requires  the  integrals  of  /  over  all  lines,  not 
only  those  through  x  or  close  to  x. 

Over  the  years  many  publications  have  dealt  with  (i),  for  example  see  [4], 
[10],  or  [22]  for  detailed  expositions  and  extensive  bibliographies.  Relatively 
recently  several  publications  have  addressed  (ii),  see  [6,  7,  30,  12].  The  basic 
idea  is  to  use  the  integrals  of  /  over  lines  close  to  x  not  to  recover  the  value 
f{x)  but  to  obtain  certain  high  frequency  content  of  /  near  x. 

A  currently  fashionable  tool  for  analyzing  the  frequency  content  of  signals 
or  functions  is  the  so-called  wavelet  transform.  For  a  sufficiently  well  behaved 
scalar  valued  function  /  defined  on  the  plane,  IR^,  its  wavelet  transform  is  a 
one  parameter  family  W^f{a,  x),  a  >  0,  of  functions  defined  by  the  formula 

(2)  W^f{a,  x)  =  \f  ^{y/a)f{x  -  y)dy 

a  JR^ 

where  ^  is  an  integrable,  usually  radial,  function  with  mean  value  0.  The 
function  is  sometimes  referred  to  as  the  analyzing  wavelet^.  The  effec¬ 
tiveness  of  the  wavelet  transform  in  practical  applications  is  well  known,  for 
example  see  [3]. 

Since  any  nice  radial  function  ^  enjoys  a  convenient  representation  as  a 
uniform  sum  of  ridge  functions,  see  [20],  the  corresponding  wavelet  transform 
of  /  can  be  computed  directly  from  the  values  of  the  Radon  transform  of  / 
via  a  formula  involving  convolution  and  backprojection^.  Such  formulas, 
unlike  (1),  involve  proper  integrals  nevertheless,  like  (1),  these  formulas  are 
generally  not  local.  One  of  the  objectives  of  this  article  is  to  show  how  such 
formulas  can  be  localized. 

^The  Radon  transform  and  its  inversion  is  discussed  in  more  detail  in  Subsection  2.1. 
A  reproduction  of  Radon’s  original  paper  [25]  is  in  [8]  and  its  English  translation  is  in  [4]. 
A  contemporary  derivation  of  this  formula  together  with  an  outline  of  Radon’s  proof  of 
its  validity  can  be  found  in  [22]. 

^Note  that  the  transform  depends  on  the  analyzing  wavelet  A  popular  choice  is  the 
Laplacian  of  the  Gaussian  or  the  so-called  “Mexican  hat” 

4'(a;)  =  (jx)^  - 

but,  of  course,  many  other  choices  axe  possible.  The  normalizing  factor  preceding  the 
integral  may  involve  a  different  power  of  a  depending  on  the  usage  and  the  authors’  tastes. 
This  definition  and  some  of  its  consequences  are  discussed  in  more  detail  in  Subsection  2.5 

®This  is  explained  in  more  detail  in  Section  2.5 


2 


1.2  Overview 


We  consider  the  natural  approximations 


(3) 


1 

gj{x)  =  --  / 

TT  Je 


dF^{q) 

q 


of  (1)  parametrized  by  positive  e  and  observe  that 
(4)  Qefix)  =  /  k{y/e)f{x  -  y)dy 


where  A:  is  a  positive  radial  function  with 


=  1  . 


Identity  (4)  means  that  gtf{x)  is,  in  effect,  the  convolution  of  /  with  an 
approximation  of  the  identity  kernel  of  “thickness”  e.  This  routinely  im¬ 
plies  that  lime_o  gtf{x)  =  f{x)  for  suitably  well  behaved  /  and,  incidentally, 
provides  an  alternate  proof  of  (1). 

Formulas  like^  (3)  and  (4)  are  the  basis  of  our  “transforms  for  quasilocal 
tomography”.  Namely,  ^(a:)  =  2k{2x)  —  k{x)  has  mean  value  0  and  the 
corresponding  wavelet  transform  can  be  computed  via 


(5) 


Wis,f{a,x)  =  f 

TT  Ja/2 


dFa,{q) 

q 


In  other  words,  the  value  of  the  wavelet  transform  Wqff{a,x)  can  be  com¬ 
puted  using  only  the  integrals  of  /  over  lines  which  intersect  the  disc  of  radius 
a.  Furthermore,  for  any  positive  e, 


(6)  E  W^f{2U,x)  =  f(,x) 

j=-oo 

and,  since  for  sufficiently  large  j  the  terms  W /(2-^e,  x)  contain  low  frequency 
information  which  can  often  be  ignored,  formula  (6)  can  be  used  as  the  basis 
of  a  “quasilocal”  approximate  inversion  formula. 

The  inversion  method  suggested  by  (5)  and  (6)  is  potentially  of  some 
practical  significance.  For  instance,  we  show  how,  in  one  particular  realiza¬ 
tion,  (5)  can  be  approximately  evaluated  using  parallel  beam  X-ray  data  and 
a  slightly  modified  version  of  the  classical  backprojection  algorithm  intro¬ 
duced  by  Shepp  and  Logan  [28]. 

In  our  development  we  include  the  following: 

•  a  brief  review  of  some  elementary  facts  concerning  the  Radon  transform 
and  its  inversion; 

^We  use  slightly  different  variants  of  the  inversion  formula  in  our  development. 
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•  a  recollection  of  the  notions  of  “local  tomography”  and  “pseudolocal  to¬ 
mography”  as  found  in  [6,  7,  30]  and  [12]  respectively  and  an  indication 
of  a  connection  between  the  two  notions; 

•  a  derivation  of  relationships  like  those  implied  by  the  pairs  of  equations 
(3),  (4)  and  (5),  (6)  and  an  indication  of  some  of  the  consequences  of 
these  formulas; 

•  a  brief  review  of  the  definition  and  the  motivation  of  the  wavelet  trans¬ 
forms,  a  general  paradigm  for  computing  such  transforms  in  terms  of 
Radon  transform  data,  and  an  indication  of  the  fact  that,  in  effect, 
“local”  and  “pseudolocal  tomography”  involve  the  computation  of  a 
wavelet  transform; 

•  discrete  analogues  of  formulas  like  (5)  and  (6)  which  lead  to  various 
reconstruction  algorithms; 

•  numerical  examples  illustrating  potential  applications  of  these  algo¬ 
rithms. 

1.3  Remarks,  conventions,  and  acknowledgements 

The  Radon  transform  and  the  notions  of  wavelets  or  wavelet  transform  have 
been  recently  treated  together  in  many  articles,  for  example  see  [1,  5,  11, 
26].  However,  to  my  knowledge,  this  is  the  first  time  that  certain  basic 
relationships  such  as  those  expressed  by  the  pairs  of  equations  (3),  (4)  and 
(5),  (6)  have  been  recorded. 

Most  of  the  development  can  be  found  in  Section  2  where,  to  maintain 
readability,  certain  technical  details  and  other  pertinent,  but  not  essential, 
remarks  are  kept  to  a  minimum.  The  main  development  is  contained  in 
Subsections  2.1,  2.3,  2.5,  and  2.7.  Subsections  2.2,  2.4,  and  2.6  contain 
material  on  “local”  and  “pseudolocal  tomography”  and  are  not  critical  to  this 
development;  they  are  included  to  bring  attention  to  the  natural  connection 
between  these  notions  and  wavelet  transforms.  Discrete  analogues  of  some 
of  the  formulas  found  in  2.3  and  2.7  are  indicated  in  Subsection  2.8.  Such 
analogues  can  be  used  as  the  basis  for  algorithms  for  the  reconstruction 
of  f(x)  and  its  wavelet  transform  from  its  Radon  transform  data;  this  is 
illustrated  by  the  numerical  examples  in  Subsection  2.9. 

Subsection  3.n  is  devoted  to  the  details  and  remarks  omitted  in  Subsec¬ 
tion  2.n. 

We  use  standard  terminology,  notation,  and  conventions.  Here  we  simply 
remind  the  reader  of  the  following: 
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•  The  convolution  f*goi  two  scalar  valued  functions  /  and  g  on  IR^  (or 
M)  is  defined  by 

/  *  9(x)  =  j  -  y)dy 

where  the  integral  is  taken  over  all  of  IR^  (or  IR)  whenever  the  integral 
is  well  defined  and  distributionally  otherwise. 

•  The  Fourier  transform  /  of  a  function  /  on  IR?  (or  IR)  is  defined  by 

KO  =  j 

where  the  integral  is  taken  over  all  of  IR^  (or  IR)  whenever  the  integral 
is  well  defined  and  distributionally  otherwise.  (^,  x)  denotes  the  usual 
scalar  product  of  ^  and  x. 

•  Whether  the  convolution  or  Fourier  transform  is  to  be  interpreted  in 
the  bivariate  or  univariate  sense  should  be  clear  from  the  context. 

•  Generic  constants,  whose  meaning  should  be  clear  from  the  context, 
are  denoted  by  c. 

In  what  follows  /  always  denotes  an  integrable  scalar  valued  function  on  the 
plane  IR^,  that  is,  fp2  \f{x)\dx  is  finite.  Other  restrictions  on  /  will  be  given 
as  needed. 

This  work  was  partially  motivated  by  the  appearance  of  [12]  which  was 
brought  to  my  attention  by  Dr.  Arje  Nachman  of  the  AFOSR. 

2  Development 

2.1  The  Radon  transform  and  an  inverse 

Recall  that  the  Radon  transform  Rf{6,  t),  0  <  6  <  27r,  — oo  <  t  <  oo,  of  an 
integrable  function  /  on  may  be  defined  by 

/OO 

f{tu0  +  sv0)ds 

-OO 

where  uo  =  {cos 6, sin 9)  and  V0  =  {-sin 6, cos 9).  Since  it  is  often  convenient 
to  view  Rf{9,  t)  as  a  family  of  functions  of  t  parametrized  by  9  we  use  the 
abbreviated  notation 

f0{t)  =  Rf{9,t)  . 

The  precise  version  of  Radon’s  inversion  formula  is 
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where 


1  /■2’f 

^x{t)  =  7r  ’ 

Ztt  Jo 


see  [25,  formula  III'].  Radon  showed  the  validity  of  this  formula  under  the 
assumption  that  /  is  continuous  and  satisfies  two  other  technical  conditions. 
Under  these  conditions  integration  by  parts  shows  the  equivalence  of  (1)  and 

(7). 

In  Subsection  2.3  we  show  that  (7)  remains  valid  under  more  general 
conditions. 


2.2  Local  tomography:  low  and  high  band  pass  filters 

Note  that  in  order  to  determine  f{x)  via  (7)  knowledge  of  /^(t)  for  all  t  and 
6  is  required.  On  the  other  hand  applying  the  formal  adjoint  of  R  to 
fg{t)  results  in 

f27V 

(8)  R*Rf{x)=  fe{{x,ue))de . 

J  0 

Observe  that  the  above  formula  for  R*Rf{x)  in  terms  of  fe{t)  is  very  local 
in  the  sense  that  only  integrals  of  /  over  lines  through  x  are  used  to  compute 
R*Rf{x).  The  relationship  between  R*Rf  and  /  may  be  more  transparent 
from 

R*Rf{x)  =  ATTJ*f(x) 


where 


J{x)  = 


Since  the  Fourier  transform  of  J*/  is 

j~m  =  i{r'/«) . 

R*Rf  may  be  regarded  as  a  low  band  pass  filtered  version  of  /. 

If  /  is  sufficiently  well  behaved  and  the  Laplacian  A  is  applied  to  both 
sides  of  (8)  then  interchanging  the  order  of  integration  and  differentiation  on 
the  right  hand  side  results  in 


r2'K 

I^R*Rf{x)=  /  fj({x,ug))de  , 

Jo 


where  fg{t)  is  the  second  derivative  of  fe{t)  with  respect  to  t.  Note  that 
formula  (9)  for  AR*Rf{x)  in  terms  of  fe{t)  is  only  slightly  less  local  than 
(8)  in  the  sense  that  only  integrals  of  /  over  all  lines  which  pass  through  an 
arbitrarily  small  neighborhood  of  x  are  needed  to  compute  AR*Rf{x). 

Using  A  to  denote  the  transformation  defined  by 

Af{x)  =  -AJ*f{x) 
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we  see  that  it  may  also  be  expressed  as 


Af{x)  =  -^AR*Rfix) 

47r 

and  computed  via  (9).  Note  that  mapping  /  — »■  J*/  is  the  inverse  of  A  since 
J^cA/  =  /,  in  other  words 


A  ^f{x)  =  J*f{x). 

Since  the  Fourier  transform  of  A/  is 

a7({)  =  \(\m . 

A/  may  be  regarded  as  a  high  band  pass  filtered  version  of  /.  Thus  one 
may  obtain  both  low  and  high  band  pass  versions  of  /,  A~^f  =  ^R*Rf 
and  A/  =  ^AR*Rf,  from  local  Radon  transform  data  via  formulas  (8) 
and  (9).  Both  versions  can  be  used  to  obtain  information  on  /  from  its 
Radon  transform  data.  For  example,  in  [6,  7,  30]  the  authors  suggest  the 
reconstruction  of  linear  combinations  of  A/  and  A~^f  rather  than  the  density 
function  /  itself;  they  refer  to  this  as  “lambda”  or  “local  tomography” . 


2.3  Low  band  pass  filters  and  inversion  formulas 

Consider  the  natural  approximation  g^f  of  (7)  parametrized  by  e  and  defined 

(10)  = 

The  fact  that  the  transformation  f  gef  is  translation  invariant  implies  that 
gef  is  the  convolution  of  /  with  some  distribution.  Appropriate  calculations 
show  that 

(11)  gj{x)  =  ke*  f{x) 

where 


k(x)  = 


ke{x)  =  e  ‘^k  (x/e)  , 

1 


•K^\X 


2^ 


1 


x(kl) , 


and 


m  ^  /  0  if  W  <  1 

\  1  if  lt|  >  1 . 


Note  that  k  enjoys  the  following  properties: 

i.  fe  is  a  non-negative  radial  function  on  IR^. 

ii.  k{x)dx  =  1 . 
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iii.  /jj2  \x\°‘k{x)dx  is  finite  for  0  <  a  <  1. 


iv.  |a:|A:(a;)  <  c\x\  ^  for  |a:|  >  2. 

Thus  ge  is  simply  the  convolution  of  /  with  a  radial  summability  kernel,  or 
approximation  of  the  identity,  of  “thickness”  e.  Indeed  formula  (11)  routinely 
imphes  (7)  whenever  /  is  continuous.  For  other  standard  consequences  of 
identity  (11)  see  the  Theorem  below. 

Use  of  the  representation® 

(12)  f(x)=^  rr 

StT^  Jq  J-oo 

together  with  the  approximation  suggested  by  (10)  results  in  a  similar  conclu¬ 
sion  with  a  different  summability  kernel.  To  wit,  consider  the  approximation 
Gif  of  (12)  defined  by 

f(T\  -1/^"/’  fe{{^^'^e)+t)-2fg{{x,ue))+fe{{x,ue)-t) 

then  the  analogous  manipulations  which  gave  (11)  from  (10)  result  in 

(14)  Gif{x)  =  Ki*f{x) 

where 


Ki(x)  =  e  {x/e)  , 


and  X  is  the  indicator  function  of  {t  :  |t|  >  1}  as  above.  The  kernel  K  may 
be  re-expressed  as 


K(x) 


TT^ 


{i-x(N)}  + 


x\'^  ||a;|  -1-  sj\x\^  -  l| 


x(kl) 


to  see  that  it  is  an  integrable  radial  function  which  is  decreasing  as  a  function 
of  |rr|  and  K{x)  =  0(|a;|~®)  as  |a;|  tends  to  oo. 

Note  that  K  enjoys  properties  i-iv  enjoyed  by  k  above.  Thus  G«/,  like 
gif,  is  the  convolution  of  /  with  a  radial  summability  kernel  of  thickness  e. 
Because  K  is  dominated  by  an  integrable  radially  decreasing  function,  namely 
itself,  it  is  significantly  better  than  k  which  does  not  enjoy  this  property. 
Relationship  (14)  together  with  the  properties  of  K  routinely  imply  various 
results  concerning  the  convergence  of  Gif(x)  to  f(x)  as  e  goes  to  0.  For 
example: 

®This  inversion  formula  can  be  found  in  [20]  on  page  196.  See  also  Subsection  3.3. 
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Theorem 


i  limGe/(a;)  =  f{x)  almost  everywhere. 

ii  GJIIli  =  0. 

iii  If  f  isin  1/  for  some  p,  1  <p  <  oo,  then  lim ||/  -  G«/||iP  =  0 . 

iv  If  /  is  in  L°°  and  continuous  at  x  then 

^^GJ{x)  =  f{x). 

Furthermore,  if  /  is  uniformly  continuous  then  this  holds  uniformly  in 

X. 


V  If  /  is  L°°  and  Holder  continuous®  at  x  of  order  a  for  some  positive  a  then 


(15) 


\f{x)-GJ{x)\  <c{ 


e“ 

e{l  +  |logc|} 


if  0  <  a  <  1, 

if  0  =  1, 

if  a  >  1 


where  c  is  independent  of  e.  Furthermore,  if  /  is  uniformly  Holder 
continuous  on  IR^  then  this  estimate  is  valid  for  all  x  with  a  constant 
c  independent  of  e  and  x. 


vi  Suppose  r  is  an  analytic  arc,  a;  G  F  and  B  is  an  open  neighborhood  of  x 
such  that  the  complement  of  F  fl  H  in  H  has  two  components  Bi  and  H2. 
If  the  restrictions  of  /  to  Bi  and  B2  have  extensions  to  the  closures  of 
Bi  and  B2  respectively  which  are  uniformly  Holder  continuous  of  order 
a  for  some  positive  a  then  (15)  holds  with  the  left  hand  side  of  the 
inequality  replaced  with 

\[  lim  /(?/)+  lim  f{y) 

^  and  y^Bi  and  y^B2 

®The  function  /  is  said  to  be  Holder  continuous  of  order  a  at  a:  if  the  k-th.  order 
diiference  of  /  in  y  at  a;  is  dominated  by  a  constant  times  |y|“,  where  k  is  the  least  integer 
greater  than  a.  Thus,  if  a  is  in  the  range  0  <  a  <  1,  this  means  that 

\f{x  +  y)-f{x)\  <  CIj/I"  , 

where  C  is  a  constant  independent  of  y;  if  a  is  in  the  range  1  <  o:  <  2,  this  means  that 
\fix  +  y)-  2f{x)  +  f{x  -  y)\  <  , 

etc.  The  function  /  is  said  to  be  uniformly  Holder  continuous  of  order  a  if  it  is  Holder 
continuous  of  order  a  at  every  point  x  and  the  constant  C  is  also  independent  of  x.  Note 
that  if  /  is  differentiable  then  it  is  H51der  continuous  of  order  one,  if  /  is  twice  differentiable 
then  it  is  Holder  continuous  of  order  two,  etc. 
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vii  If  /  is  in  2/  for  some  p,  2  <  p  <  oo,  then  G^f  is  Holder  continuous  of 
order  1  —  2/p.  If  /  is  Holder  continuous  of  order  ol  then  GJ  is  Holder 
continuous  of  order  a  +  1. 

Items  ii  -  vii  remain  valid  if  is  replaced  with  pe/.  See  Subsection  3.3 
for  more  details.  Similar  results  are  valid  for  regularized  variants  of  GJ  and 

gj- 

Note  that  the  Fourier  transform  of  Ge  is  given  by 

Gj({) = Kiom  ■ 

Since  K  is  continuous,  K{0)  =  1,  and  lim|^|_»c)o-f^(0  =  0,  G^  is  the  result  of 
processing  /  through  a  low  frequency  band  pass  filter  which  converges  to  the 
identity  as  e  tends  to  0.  Similar  remarks  ar^aUd  concerning  g^f  with  ^(e^) 
playing  the  role  of  the  low  band  pass  filter  K{e^). 


2.4  Local  or  pseudolocal  tomography 

For  sufficiently  well  behaved  /  the  inversion  formula  (1)  can  be  re-expressed 


as 

(16) 


=  L - i - 


where  /^(t)  is  the  derivative  of  fe{t)  with  respect  to  t.  In  view  of  this  expres¬ 
sion  the  authors  of  [12]  consider  the  decomposition 


(17) 


/(l)  =  hj(x)  +  {/(l)  -  hj{x)} 


where 

(18) 


t 


and  suggest  the  reconstruction  of  the  quantity  hef  rather  than  the  density 
function  /  itself  to  obtain  high  frequency  information  concerning  /.  They 
refer  to  this  as  “pseudolocal  tomography” . 

Note  that  (17)  may  be  expressed  as 


f{x)  =  hj{x)  +  gj{x) 


where 

(19)  hj{x)  =  f(x)  -  gj(x) 

and  g^f  is  defined  by  (10).  The  last  relation  implies  that  the  Fourier  trans¬ 
form  of  hef  may  be  expressed  as 

OiO  =  {1  -  i(«fl}/(0 
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where,  as  indicated  in  the  previous  subsection,  is  a  low  frequency  band 
pass  filter.  Thus  h^f  may  be  regarded  as  the  result  of  processing  /  through 
a  high  frequency  band  pass  filter  which  converges  to  0  as  e  tends  to  0. 

A  decomposition  analogous  to  (17)  can  be  made  using  representation  (12). 
Namely 

f{x)  =  HJ{x)  +  GJ{x) 

where 

HJ{x)  =  f(x)  -  GJ(x) 

and  Gef  is  defined  by  (13).  In  this  case,  for  sufficiently  well  behaved  /,  the 
formula  for  Hif(x)  may  be  expressed  as 

(20)  H.f(x)  =  ^ 

It  should  be  clear  that  Hef  may  also  be  regarded  as  the  result  of  processing 
/  through  a  high  frequency  band  pass  filter  which  converges  to  the  0  as  e 
tends  to  0. 

It  is  important  to  note  that,  by  virtue  of  (18)  and  (20),  both  htf{x)  and 
Hef{x)  are,  in  principle,  computable’^  in  terms  of  integrals  of  /  over  lines 
which  intersect  the  disk  of  radius  e  centered  at  x.  For  this  reason  they  may 
be  regarded  as  being  computable  from  local  Radon  transform  data  of  /. 

Thus  “local  tomography”  whose  objective  is  to  compute  A/  and  “pseu¬ 
dolocal  tomography”  whose  objective  is  to  compute  hef  appear  to  take  dif¬ 
ferent  routes  to  roughly  the  same  goal.  However,  the  fact  that  both  A/  and 
hef  are  high  frequency  band  filtered  versions  of  /  which  can  be  computed  in 
terms  of  local  Radon  transform  data  is  not  the  only  similarity  between  these 
functions.  Observe  that  the  integrand  in  (18)  may  be  replaced  by 

/g((^,^g)  -  i)  -  feii^^ue)) 
t 

which  in  turn  is  well  approximated  by  —fe{{x,U0))  if  c  is  sufiiciently  small. 
Thus  for  sufficiently  small  e  we  may  write 

hef{x)  j4  /  ((^.  '^e))dQ  =  -A/(a:) 

47r‘^  Jq  tt 

or,  more  precisely,  for  sufficiently  well  behaved  / 

(21)  hef{x)  =  —Af{x)  +  0(e)  as  e  0. 

TT 

Further  connections  are  discussed  in  Subsection  2.6. 

^These  formulas  are  valid  if  /  is  sufficiently  smooth.  Otherwise  one  should  use  “regu- 
larized”  analogues,  see  Subsection  2.6. 
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2.5  Wavelet  transforms  and  tomography 

Recall  that  the  wavelet  transform  of  /  via  the  “wavelet”  ^  is  defined  by 

Wysff{a,x)  = 

where  a  >  0,  =  a~^'if{x/a),  and  ^  is  an  integrable  function  or  measure 

on  with  mean  value  0  or,  in  other  words,  ^(0)  =  0.  Note  that  by  virtue 
of  Plancherel’s  formula 

Wii,f(a,x)  =  (27r)"2  /  • 

JR^ 


So  if  the  support  of  4'  is  the  annulus  0  <  bo  <  \^\  <  bi  <  oo,  this  is 
“roughly”  the  case  for  some  bo  and  bi  whenever  ^  is  a  smooth  and  integrable 
radial  function  with  mean  value  zero,  then  W^f{a,x)  may  be  viewed  as  the 
frequency  content  of  /  in  the  band  6o/a  <  |^|  <  fei/a  at  a;. 

In  theoretical  applications,  for  example  see  [2, 14,  24],  the  specific  ^  is  not 
really  important,  only  certain  of  its  properties  play  a  major  role.  However, 
in  numerical  work  it  is  nice  to  have  explicit  expressions  for  both  ^  and  its 
Fourier  transform.  In  addition  to  the  Mexican  hat  wavelet  mentioned  earlier, 
typical  examples  are  the  following: 

i  Derivatives  of  the  Poisson  kernel 

D  a 

^  27r(a2  +  |a:|2)3/2 

whose  Fourier  transform  is 


For  instance 


dPa  ,  ,  2a^  —  a\x\^ 

^  da  ^  27r(o2  + 


whose  Fourier  transform  is 


. 


ii  Appropriate  differences  of  known  summability  kernels  or  measures.  For 
instance 

^o(a:)  =  6(a;)  -  Pa{x) 


or 

^a(a:)  =  2Pa(2x)  -  Pa{x) 


12 


where  Paix)  is  the  Poisson  kernel  and  6(x)  =  a  ‘^d{x/a)  is  the  bivariate 
normalized  Dirac  delta  “function” .  In  this  case 

f  =  1  - 

or 

respectively. 

If  ^  is  a  radial  function  which  is  sufficiently  smooth  then  it  can  be  rep¬ 
resented  as  a  uniform  sum  of  ridge  functions,  that  is 

1 

(22)  i;{{x,U0))de 

for  some  even  univariate  function  ip(t).  The  significance  of  representation 

(22)  lies  in  the  fact  that,  for  fixed  6,  the  (bivariate)  convolution  of  tp({x,  ue)) 
with  f{x)  can  be  expressed  as  the  (univariate)  convolution  of  with  fg{t) 
evaluated  at  {x,ue).  This  leads  to 

(23)  *  fi^)  =  7^ -  /  i’a*  fe{{x,  Ug))de 

where  V’o(^)  =  In  short,  (22)  gives  rise  to  a  convenient  formula, 

namely  (23),  for  computing  the  wavelet  transform  of  /  with  wavelet  W  in 
terms  of  the  Radon  transform  of  /. 

If  ^  and  •0  satisfy  (22)  there  are  various  alternate  formulas  relating  ^  to 
0,  see  [20].  A  particularly  convenient  one  relates  their  Fourier  transforms, 

(24)  2^(|^|)  =  |e|4’(0. 

whenever  ^  is  a  smooth  function  which  decays  sufficiently  rapidly.  Among 
other  things,  this  formula  implies  that  if  ^  is  sufficiently  smooth  then  0  is 
integrable,  has  mean  value  zero,  and 

0a  *  fe{t) 

is  a  (univariate)  wavelet  transform  of  fg.  Thus  in  this  case  (23)  impfies  that 
the  (bivariate)  wavelet  transform  of  /  is  a  sum  of  appropriately  backprojected 
(univariate)  wavelet  transforms  of  the  /$’s,  in  other  words, 

(25)  W^f{a,x)  =  - — /  W^f0{a,{x,ug))dd . 

2i7r a  Jo 

For  example,  the  bivariate  wavelet 

.  2a®  -  a\x\^ 

~  27r(a2  +  |a:|2)5/2 
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mentioned  earlier  is  related  via  (22)  to  the  univariate  wavelet 


ipa{t)  = 


a®  -  3aH^ 

'K{a?‘  +  ^2)3 


with 


7)_  (  t\ 


a2^2g-|ar| 


Other  examples  are  included  in  what  follows. 


2.6  A  connection  with  local 

and  pseudolocal  tomography 

According  to  the  authors  of  [6,  7,  30],  in  practical  applications  of  local  to¬ 
mography  they  do  not  compute  Af(x),  rather  they  attempt  to  reconstruct 
A($e  *  f){x)  for  some  approximation  of  the  identity  $e(a;)  =  €~‘^^(x/e), 
e  >  0,  where  $  is  an  integrable  function  with  total  integral  one.  Now,  for 
sufficiently  well  behaved  $,  A$  is  an  integrable  function  with  mean  value 
zero,  A($e)  =  e“^(A$)e,  and  we  may  write 

(26)  A($£*/)(a;)  =  e“^(A$)e*/(a:). 

The  last  formula  may  be  re-expressed  as 

(27)  A{^,*f){x)  =  e-^W^f{e,x) 
where 

(28)  ^(a;)  =  A$(a;) . 

In  other  words,  eA($j  *  f){x)  with  a  =  e  is  a  wavelet  transform  of  /. 

If  $  is  a  sufficiently  well  behaved  radial  function  then,  in  view  of  (25),  the 
wavelet  transform  suggested  by  (27)  can  be  conveniently  computed  in  terms 
of  the  Radon  transform  of  /.  Indeed  the  corresponding  univariate  wavelet  tp 
can,  in  certain  cases,  be  easily  determined  from  $.  For  example,  by  virtue 
of  (24), 

(29)  2^(|f|)  =  |?|^»(f) 

whenever  relations  (22)  and  (28)  are  valid  and  'Jf  is  sufficiently  well  behaved. 
In  particular,  if 

^  27r(l  -f  \x\^YA 

then  $£(a;)  is  the  so-called  Poisson  kernel  for  the  upper  half  space  {(x,e)  : 
X  G  and  0  <  e  <  (X)},  (28)  gives 

~  27r(l  +  |a;|2)5/2  ’ 
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and  the  relationship  between  ^  and  '0  results  in  the  formula 

, .  ,  1-3^2 

This  is  essentially  the  example  considered  in  Subsection  2.5.  The  case 

p-klV2 

$(a;)  = 


27r 


gives  rise  to 

so  the  Gaussian  kernel  gives  rise  to 

(x,  u,))M 

zTre  Jo 

where  V'  is  the  univariate  Mexican  hat  wavelet. 

The  expression  h^f{x)  from  the  “pseudolocal  tomography”  found  in  [12] 
is  a  wavelet  transform  of  /  with  a  =  e  and 


=  S(x)  -  ke(x) 

where  ke(x)  is  described  by  the  formulas  immediately  succeeding  (11)  and 
S(x)  is  the  bivariate  normalized  Dirac  delta  “function”.  Thus  this  analyzing 
wavelet  is  a  special  case  of  the  general  form  mentioned  in  example  ii  in 
previous  subsection. 

On  the  other  hand,  the  authors  of  [12]  also  suggest  the  computation  of 

*/)W 

for  positive  ei  and  e2  in  lieu  of  the  computation  of  hef{x).  Since,  by  virtue 
of  (19)  and  (11), 

h«i($£2  *  f)ix)  =  ($£2  -  *  f{x) 

where  ke^{x)  is  the  kernel  in  equation  (11),  it  is  clear  that  for  suitably  well 
behaved  $  this  expression  is  some  sort  of  wavelet  transform.  Indeed,  if 
ei  =  ce2  =  ce  then 

heii^e2  *  f){^)  =  e~^Wii,f{e,x) 

where  ^(a:)  =  ^{x)  —  kc*^{x).  If  $  is  radial  then  this  wavelet  transform  may, 
in  theory,  be  computed  in  terms  of  the  Radon  transform  of  /  in  a  manner 
similar  to  that  suggested  above.  Furthermore,  if  e  is  sufficiently  small  then 
one  may  use  the  approximation  suggested  by  (21) 

h.,($e,*/)(a:)«-A($,*/)(a:) 

TT 

together  with  the  formulas  for  computing  A($£  * 

Of  course,  similar  remarks  are  also  valid  for  Hef{x). 
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2.7  A  wavelet  transform  for  quasilocal  tomography 


Suppose 

(30)  =  Kxj^ix)  -  Ki{x) 

where  Kg,  e  =  1/2  or  1,  is  the  convolution  kernel  in  Formula  (14).  This  ^ 


can  also  be  described 

via 

l2a;P  -  1 

1^1 

X(|2x|)  + 


\/kP-l 


rr 


x(kl) 


In  view  of  (13)  and  (14)  it  follows  that 


(31) 


^  f(x)  f  /g((rg,ug)  +t)-2fe{{x,U0))+fe{{x,ue)-t) 

“  StT^  Jo  Ja/2<\t\<a 

Prom  (31)  it  is  clear  that  the  wavelet  transform 


W<5,f{a,x)  =  <ifa*f{x) 

does  not  require  knowledge  of  the  fuU  Radon  transform  of  /,  it  can  be  com¬ 
puted  in  in  terms  of  only  those  integrals  of  /  over  lines  whose  intersection 
with  the  disk  {x  :  |3;|  <  a}  is  non-empty. 

Viewed  in  terms  of  the  material  in  Subsection  2.5  the  function  ■0  corre¬ 
sponding  to  ^  is  the  distribution 


X{2t)  -  x(t)  \ 
2t2  j 


where  S{t)  is  the  univariate  normalized  Dirac  delta  “function”  and  x(t)  is  the 
indicator  function  of  the  set  {t :  |t|  >  1}.  This  is  an  immediate  consequence 
of  (31). 

Observe  that 


N 

51  *  f{x)  =  KaMg  f{x)  -  K2N,  *  f{x) 

J=M+1 

for  any  integers  M  <  N  and  any  positive  e.  Since, 

lim  K2Mg  *  f{x)  =  f{x)  and  lim  K2n^  *  f{x)  =  0 


00 

f{x)=  ^2i€*/(a:). 

j“-0O 


it  follows  that 
(32) 
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In  view  of  (31),  formula  (32)  is  an  inversion  formula  for  the  Radon  transform. 
Truncating  (32)  appropriately  leads  to  approximate  reconstruction  formulas. 

For  example,  If  /  is  supported  in  the  unit  disk  {a;  :  |a;|  <  1}  then  the 
integrand  in  formula  (31)  reduces  to  -2fe{{x,U0))/t^  whenever  la;|  <  1  and 
\t\  >  2.  Hence  if  |a;|  <  1  it  follows  that 

1 

*  fix)  =  feiix,  U0))de 

whenever  a  is  greater  than  4.  In  this  case  formula  (32)  reduces  to 

1  ^ 

(33)  fix)  =  ^^22No- feiix, U0))de  +  ^2.'^  *  fix) 

where  Nq  is  the  least  integer  so  that  2^“  >  4/e.  In  numerical  applications 
considerations  of  resolution  impose  natural  truncation  points  for  the  last 
sum.  See  the  next  subsection. 

Remark  1  It  is  important  to  note  that  the  wavelet  ^  defined  by  (30) 
can  be  replaced  with 

K,,ix)  -  K,,ix) 

for  any  ei  and  €2  which  satisfy  0  <  Ci  <  €2  with  consequences  analogous  to 
those  indicated  above  for  the  special  case  ei  =  1/2  and  ^2  =  1.  This  flexibility 
may  be  signiflcant  in  numerical  experiments  and  practical  applications. 

Remark  2  Of  course  many  of  the  above  observations  are  also  valid  if  K 
is  replaced  by  the  kernel  k  in  formula  (11)  or  certain  other  kernels,  including 
certain  regularizations  of  k  and  K.  However  k  is  not  as  well  behaved  as  K 
and,  although  other  kernels  may  be  better  behaved,  it  is  unlikely  that  the 
formula  for  the  corresponding  wavelet  transform  will  be  as  simple  as  (31). 


2.8  Discrete  analogues 

Suppose  /  is  supported  in  the  unit  disk  {x  :  |a;|  <  1},  we  are  given  the 
discrete  Radon  transform  data 

(34)  Rfi^,^),  m  =  0,l,...,M-l,  n  =  0,±l,±2,...,±Ar, 


and  we  wish  to  approximate 


(35)  f{x) 


-1 

Stt^  Jo  7-00 


or 

(36) 


1 

w^fia,  x)  =  - —  /  'ifa*  feiix,  U0))de 
Z'Ka  Jo 
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where 


x(2t)-xW1 

j  ■ 


Here  6(t)  is  the  univariate  normalized  Dirac  delta  “function”  and  x(i)  is  the 
indicator  function  of  the  set  {t :  |t|  >  1}. 

View  the  inner  integral  in  either  formula  as  the  (univariate)  convolution  of 
an  appropriate  distribution  w  with  fe  and  approximate  the  outer  integration 
via  the  periodic  trapezoid  rule.  This  results  in  the  approximation 


27r 

(37)  c—  w  *  fe^iix,  uej) 

m=0 

where  6m  =  m  =  0, —  1,  and  c  is  the  appropriate  constant. 

There  are  several  ways  to  approximate  the  convolutions  w  *  .  We  use 

the  paradigm  suggested  in  [16].  Namely,  first  replace  w  with  the  distribution 

OO 

(38)  w{t)  =  53  'WnS{t  -  §) 


where  6{t)  is  the  univariate  unit  Dirac  delta  “function”  at  the  origin, 


/OO 

w{t)(j)n{t)dt , 

“OO 


and  (f>nit)  is  sufficiently  smooth,  is  centered  and  concentrated  around  n/N, 
and  has  total  integral  one*.  Next  replace  fg^  with 

(39)  /„(()=  E 

n=-N 

where  A  is  a  continuous  function  which  satisfies 


1  if  n  =  1 
0  otherwise 


so  that  fm{t)  is  a  continuous  function  which  interpolates®  fd^{t)  att  —  n/N, 
n  =  0,  ±1,  ±2, . . ..  Observe  that 


(40) 


W  *  fm{t) 


<x>  (  N 

13  {  tl 

n=-oo 


^Wn  is  the  distribution  w  evaluated  at  which,  instead  of  the  integral,  may  be  more 
properly  expressed  as  Wn  = 

should  also  be  chosen  so  that  fm{i)  approximates  well. 
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noting  that 


N 

l=-N 

and  if  A  has  compact  support  then  the  sum  in  n  in  (40)  has  a  finite  number 
of  non-zero  terms. 

Since  fe^{^)  =  Rf{^,§),  the  approximations  suggested  by  (37)  and 

(40)  give  rise  to  approximations  of  (35)  and  (36)  which  can  be  computed  in 
terms  of  the  discrete  Radon  transform  data  (34)  of  /. 

For  example,  if 

(41)  (f>r^{t)  =  N(j){Nt-n) 
where 

(«)  =  {  0  Otis 

then  the  distribution  w  in  (35)  gives  rise  to 

yj  -n)-  2(f>{n)  +  (pjNt  -I-  n) 

”  7-00 


which  simplifies  to 
(43) 


Wn  = 


8N^ 

—  1 


The  resulting  approximation  Af  of  /  is  given  by 


(44)  Af{x) 

In  the  case 

(45) 


AixM 


M-1  oo  ( 

N 

E  E 

l=~N  ^ 

m=0  n——oo 

m  =  I  ^ 

0 

ujmiA 

—  4(n  —  ly  j 

if  \t\  <  j! 

otherwise 


K{x,uej  -  §) 


the  algorithm  suggested  by  (44)  is  the  well  known  algorithm  introduced  in 
[28]. 

To  approximate  the  wavelet  transform  (36)  observe  that  the  distribution 
w  —  Tpa  iJi  (36)  and  the  (j)  in  (41)  give  rise  to 


Wr, 


N 


L 


(l>{Nt  —  n)  —  2(j){n)  -f  4>{Nt  +  n) 

,/2<|t|<o 


dt . 


To  simplify  evaluation  of  this  expression  we  choose  a  convenient  relationship 
between  the  parameters  a  and  A/',  namely, 


(46) 
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for  some  integer  ^  =  0, 1, 2, . . ..  If  (46)  holds  then 


(47) 


Wn 


(  4^2 


I  4n*— 1 
J  4N^ 

I 

2<+i(2<+i-l) 

I  0 


if  n  =  0 

if  <  \n\  <  2^ 
if  ]n|  =  2^-1 
if  |n|  =  2^ 
otherwise. 


Denote  the  train  of  6  functions  with  weights  given  by  (47)  as  Wa  to  indicate 

its  dependence  on  a.  Here  a  takes  on  the  possible  values  2^/N,  ^  =  0, 1, 2, _ 

Then 

_i  M-l 

(48)  AWiSff{a,x)  =  J]  *  fm{{x,UgJ)) 

m=0 

is  a  natural  approximation  to  a  wavelet  transform  of  /  which  is  computable 
in  terms  of  its  Radon  transform  data  (34).  Furthermore,  if  X{t)  is  chosen  to 
have  small  support  such  as  (45)  then  AWii,f(a,  x)  depends  only  on  integrals  of 
/  over  lines  which  intersect  the  disc  of  radius  a  centered  at  x\  more  precisely, 
AWi^f{a,x)  depends  only  on  Rf{mTT/M,n/N)  for  -  {x,Urmr/M)\  <  a  and 
m  =  0, 1, . . . ,  M  —  1. 

Finally,  if  we  denote  the  train  of  6  functions  with  weights  given  by  (43) 
as  w,  it  is  not  difficult  to  see  that 

OO 

(49)  w(t)  =  Y,W2t/N{t)  ■ 

e=o 

Relationship  (49)  is  essentially  a  discrete  analogue  of  (32).  Using  reasoning 
similar  to  that  which  led  to  (33)  results  in  the  approximation 

No  7y  M-l 

(50)  Af(x)  =  Y:  AW^nVjN,  x)  +  Y  “»))  ■ 

j=0  m=0 

where  Nq  is  the  smallest  integer  such  that  2^“  >  2N.  What  is  remarkable 
about  this  simple  decomposition  is  the  fact  that  often  the  first  few  terms 
of  the  first  sum  are  sufficient  to  obtain  a  reasonable  reconstruction.  This  is 
illustrated  by  Example  1  below. 

Remark  1  Similar  formulas  are  valid  involving  the  more  general  'F’s 
mentioned  in  Remark  1  at  the  end  of  the  previous  subsection.  The  signif¬ 
icance  of  this  lies  in  the  fact  that  one  may  want  to  experiment  with  the 
various  choices  of  parameters  ei  and  62  to  determine  selections  most  suitable 
for  specific  applications. 

Remark  2  The  choice  of  0  suggested  by  (42)  was  motivated  by  the 
connection  with  the  classical  algorithm  of  Shepp  and  Logan  [28].  However, 
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in  view  of  the  distribution  w  in  (35)  and  (36),  this  choice  seems  somewhat 
less  than  ideal. 

A  smooth,  yet  extremely  simple,  selection  is 


(51) 


if  |t|  <  1 
otherwise 


which  is  twice  continuously  differentiable,  is  supported  in  {t  :  |t|  <  1}  and 
satisfies 


OO 

0(0)  =  !>  0(^  -n)  =  l, 

n=— OO 


and 


=  1. 


In  this  case  the  distribution  w  corresponding  to  the  w  in  (35)  has  the  weights 
(52) 


Wn  =  C< 


1 

l-21og2 

l-(|n|  +  n^)log(H±l)-(|»| 


’*^)log(r*) 


if  n  =  0 
if  |n|  =  1 
if  \n\  >  2 


where  c  is  an  appropriately  chosen  constant.  In  other  words,  (52)  is  the 
analogue  of  (43)  if  (51)  is  used  instead  of  (42). 

The  weights  for  the  distributions  Wa,  a  =  2^/N,  £  =  0,1,2 which  may 
be  used  to  approximate  the  wavelet  transform  (36)  are  the  following;  when 
£  =  0  then 


Wn 


3 

-3/2 

0 


if  72  =  0 
if  |n|  =  1 
if  jnj  >  2  ; 


when  £  =  1  then 


2  if  n  =  0 

-16 -24  log  2  if|n|  =  l 

241og2-17  ifini  =  2 

0  if  jrzj  >  3 ; 


and  when  £>2  then 
(53) 


Wn  =  2c< 


1/2^ 

6{l-(|n|+n^)log(H±l)-(|n|-n^)log{pHl)} 
-l/|n|+3  +  6W-6(H+»^)log(H±l) 
l/M+3-6H-6(W-».^)log(JfL) 


0 


if  n  =  0 

if  2^-^  <  |n|  <  2^ 
if  |n|  =  2^~^ 
if  |n|  =  2^ 
otherwise. 
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where  c  is  an  appropriately  chosen  constant  and  is  the  same  in  all  three 
formulas.  In  other  words,  (53)  and  the  two  preceding  formulas^®  are  the 
analogue  of  (47)  if  (51)  is  used  instead  of  (42).  The  wavelet  transforms  in 
our  numerical  examples  are  computed  with  these  weights. 

2.9  Numerical  examples 

We  present  the  results  of  two  numerical  experiments  based  on  the  algorithms 
suggested  by  the  formulas  in  Subsection  2.8.  As  is  customary,  we  used 
the  celebrated  phantom  $  which  was  originally  introduced  and  documented 
in  [28].  That  is,  $(a:)  is  the  linear  combination  of  indicator  functions  of 
ellipses  as  described  in  [28].  MATLAB  software  was  used  to  implement  the 
algorithms  and  display  the  results. 

2.9.1  Example  1 

The  point  of  this  example  is  to  compare  a  classical  reconstruction  with  that 
suggested  by  the  discrete  analogues  of  the  wavelet  transform. 

We  used  the  data  described  by  (34)  with  f{x)  =  ^{x),  M  —  90,  and 
A”  =  100.  The  results  are  summarized  by  Figures  1-5.  Each  “reconstruction” 
was  evaluated  on  a  201  x  201  grid  determined  by  a;  =  ('^Yoo°°> 
represented  by  a  matrix  /  =  {Ij2,ji)  where  ji  and  j2  take  on  the  values 
0,1,...,  200.  Each  figure  contains  three  plots. 

•  A  plot  of  the  “cross-section”  /loiji,  ji  =  0, 1, . . . ,  200. 

•  A  plot  of  /tioiji,  ji  =  0,1,..., 200,  where  It  is  the  range  truncated 
variant  of  I  defined  by 


r  fj,i  if  Ijiji  > 

^^j2,h  ~  j  ^  ^jiJi  —  Ml 

[  1^0  if  Ih,h  <  /^o  • 

The  values  of  iiq  and  were  chosen  so  that  the  variations  inside  the 
“skull”  would  be  visible  and  were  determined  experimentally  by  viewing 
the  plot  of  hoiji  with  different  scalings  on  the  vertical  axis. 

•  A  gray  level  representation  (image)  of  Jt  using  256  gray  levels  uniformly 
distributed  between  fio  and  yUi. 

Figure  1  gives  the  results  of  the  computation  of  Af{x),  the  approximation 
of  /  given  by  formula  (44).  In  other  words,  this  is  the  result  of  applying 

Since  all  three  formulas  are  special  cases  of  one  general  formula  we  will  refer  to  them 
collectively  as  (53)  in  what  follows. 
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the  celebrated  Shepp-Logan  algorithm  [28].  It  is  presented  for  comparison 
purposes. 

We  computed  the  approximations  AW^f{2^/N,  x),  j  =  0, 1, . . . ,  7,  to  the 
wavelet  transform  using  the  algorithm  suggested  by  (48)  and  the  weights 
given  by  (53). 

Figures  2,  3,  and  4  show  AW^f{a,x)  with  a  =  1/N,2/N,  and  4/N  re¬ 
spectively.  Figure  5  shows  the  reconstruction 

7  T\T  M-1  _ 

Af(x)  =  Y,  AW^J(2l/N,  x)  +  ^  £  /„({!,  u,)) . 

m=0 

The  algorithms  for  computing  Af{x)  were  not  normalized  so  that  the 
plots  in  Figures  1  and  5  should  only  be  compared  relatively  not  absolutely. 

2.9.2  Exeimple  2 

The  point  of  this  example  is  to  numerically  illustrate  the  local  nature  of  the 
wavelet  transform  W^f{a,x)  with  ^  as  described  by  (30),  when  computed 
from  Radon  transform  data. 

We  used  the  data 

(54)  R/(^,f),  m  =  0,l,...,M-l,  n  =  0,±l,±2,...,±iV, 

with  f(x)  —  ^{x  +  xq),  M  =  90,  and  N  =  200.  Here  xq  —  (0,  -.605)  is  the 
center  of  the  central  small  ellipse  in  the  lower  portion  of  the  original  phantom 
$.  We  computed  “reconstructions”  from  from  various  truncations  of  the 
data  matrix,  Dm,n  =  second  variable;  more  specifically 

n  =  0,  ±1,  ±2, . . . ,  ±rN,  with  r  =  1, 3/4, 1/2, 1/4,  and  1/8. 

The  results  are  summarized  in  Figures  6-10  and  11-15. 

Each  “reconstruction”  was  evaluated  on  a  401  x  401  grid  determined  by 
^  represented  by  a  matrix  I  =  (/,-2,ii)  where  ji 

and  j2  take  on  the  values  0, 1, . . .  ,400.  Each  of  the  Figures  6-10  contains 
three  plots. 

•  A  gray  level  representation  (image)  of  the  data  using  256  gray  levels 
uniformly  distributed  between  0  and  the  maximum  of  the  data. 

•  A  gray  level  representation  (image)  of  the  truncated  matrix  It,  with 
values  between  and  which  represents  the  result  of  simply  apply¬ 
ing  the  Shepp-Logan  algorithm  (formula  (44))  to  the  data;  256  gray 
levels  uniformly  distributed  between  yUo  and  were  used.  The  range 
truncated  matrix  It  was  obtained  by  essentially  the  same  procedure 
as  described  in  Example  1;  namely  the  values  of  //q  and  /ii  were  cho¬ 
sen  so  that  the  variations  inside  the  “skull”  would  be  visible  and  were 
determined  experimentally  by  viewing  the  plot  of  /2oi,ji  with  different 
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scalings  on  the  vertical  axis.  See  Figures  11-15.  In  the  case  r  =  1/8  it 
was  impossible  to  find  appropriate  values  of  Ho  and  fii,  the  case  r  =  1/4 
was  at  best  borderline. 

•  A  gray  level  representation  (image)  of  the  truncated  matrix  It,  with 
values  between  //q  and  /xi,  which  represents  AW^/(l/Ar,  x)  as  computed 
by  the  algorithm  suggested  by  (48)  and  the  weights  given  by  (53);  256 
gray  levels  uniformly  distributed  between  jjio  and  were  used.  The 
range  truncated  matrix  It  was  obtained  as  described  in  the  previous 
item.  Note  that  even  in  the  cases  r  =  1/4  and  1/8  it  was  possible  to 
find  appropriate  values  of  /xq  and  fii. 

Figures  6-8  are  scaled  so  that  they  occupy  the  same  area.  Figures  8-10  are 
scaled  in  correct  proportion. 

Figures  11-15  contain  the  central  “cross-sections”  of  the  reconstructions 
found  in  Figures  6-10  respectively  and  document  how  the  range  truncated 
matricies  It  were  obtained.  Each  figure,  call  it  by  it’s  number  figure  number, 
consists  of  four  plots  arranged  as 


1  3 

2  4 

and  which  contain,  respectively,  the  following. 

1.  A  plot  of  the  “cross-section”  hoiji,  ji  =  0, 1, . . . ,  400,  where  the  401  x 
401  matrix  I  represents  the  result  of  simply  applying  the  Shepp-Logan 
algorithm  (formula  (44))  to  the  corresponding  data  in  Figure  n  —  5. 

2.  A  plot  of  7^201, ii)  ji  =  0,1,..., 400,  where  It  is  the  range  truncated 
variant  of  the  matrix  I  defined  in  item  1.  As  mentioned  earlier,  this 
truncation  was  obtained  in  the  same  manner  as  indicated  in  Example 
1. 

3.  A  plot  of  the  “cross-section”  hoiji,  ji  =  0,1, . . . ,  400,  where  the  401  x 
401  matrix  /  represents  AW^f{l/N,x)  as  computed  from  the  data  in 
Figure  n  —  5  by  the  algorithm  suggested  by  (48)  and  the  weights  given 
by  (53). 

4.  A  plot  of  7^201  ji)  ji  —  0,1,..., 400,  where  It  is  the  range  truncated 
variant  of  the  matrix  7  defined  in  item  3.  Again,  this  truncation  was 
obtained  in  the  same  manner  as  indicated  in  Example  1. 
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3  Details 


3.1 


To  see  that  Fx{t)  is  the  average  of  the  integrals  of  /  over  lines  which  are  a 
distance  |f|  from  x  simply  write 


rZTT  rcx> 

2TrF^{t)  =  /  f{{{x,  ue)  +  t]u0  +  sv0)dsd6 

jO  J — oo 
r27r  roo 

=  /  f{{{x,U0) +  t}ue  + {{x,vo)  +  s}v0)dsd0 

Jo  J — oo 


because  the  integration  in  s  is  translation  invariant.  Since  x  =  {x,  ue)u0  + 
{x,V0)v0,  we  may  write 


which  is  the  desired  result  since  the  expression  in  braces  is  the  integral  of  / 
over  the  line  which  is  t  units  away  from  x  in  the  direction  U0. 

3.2 

The  formal  adjoint  of  R  is  the  transformation  which  maps  suitably  well 
behaved  scalar  valued  functions  f*{9,  t)  on  [0, 27r)  x  ]R  to  functions  R*f*(x) 
on  and  satisfies  the  relation 

j  j{x)R*nx)dx=  r  r  Rf{e,t)r{e,t)dtde. 

JR^  JQ  J-oo 

Routine  calculations,  namely 

f2TV  POO  p2'K  poo  poo 

/  /  Rf(6,t)f*{6,t)dtdd=  /  /  f{tu0  + sv0)f*{O,t)dsdtdO 

Jo  J —oo  Jo  J-ooJ—oo 


p2n  p 

=  fix)f*{0,  {x,  U0))dxd9 

JQ  JR? 

where  t  =  (re,  uo)  and  s  =  (re,ne),  show  that 

R*f*(x)=  [  f*{e,{x,U0))de . 

Jo 

To  see  that  R'^Rf  =  47rJ*/  write 

p2'k  p2'K  poo 

R*Rf{x)=  fe({x,U0))d6=  /  f{{x,U0)u0  +  sv0)dsd6 . 

Jo  Jo  J —oo 
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Because  the  inner  integral  is  translation  invariant  the  variable  s  may  be 
replace  by  {x,  vg)  +  s  which,  in  view  of  the  fact  that  x  =  {x,  U0)ue  +  (x,  ve)v0, 
results  in 

/•27r  POO 

R*Rf{x)  =  /  f{x  +  sve)dsd9 

Jo  J  ~oo 

=  2  r  r 

Jo  Jo  5 

= L 

when  the  polar  change  of  variables  y  =  sue  is  used. 

The  fact  that  the  Fourier  transform  of  |a:|"^  is  a  constant  multiple  of 
itself  follows  from  elementary  fact  that  the  bivariate  Fourier  transform  of  a 
radial  distribution  homogeneous  of  degree  minus  one  is  a  radial  distribution 
homogeneous  of  degree  minus  one.  The  value  of  the  constant  follows  from 
the  Fourier  inversion  formula 

Recall  that  the  Laplacian  A  is  the  differential  operator  defined  by 
A/(x)  =  +  ||(x) 

whose  Fourier  transform  is 


A/(0  =  -l^lV(0- 

In  view  of  this  all  the  formulas  involving  A  and  are  routine  consequences 
of  previous  identities. 


3.3 

Formula  (12)  is  not  equivalent  to  (1).  However  it  is  a  natural  consequence  of 
certain  summability  formulas  for  the  inversion  of  the  Radon  transform.  The 
details  can  be  found  in  [20,  page  196]  and  [21]. 

The  formulas  gef{x)  =  ke*f{x)  and  Gef{x)  =  Ki*f{x)  follow  from  the 
identities 

/  fe{{x,U0) -t)de  =  2  f  -^j=^==dy, 

Jo  J\y\ 


(55) 

(56) 


f  f  =  -  / ,  -^Jix  -  y)dy , 

Jo  J\t\>t  e  Jr^  |y| 
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and 


(57)  /  "  / 

JQ  J\t 


f  fe{{x,  ue)  - ,  f  VbAP  -  1 


\t\>e 


=4/ 

Jh 


h>£  |yp 


f{x-y)dy. 


Let  Ii{t),  /2(e),  and  /3(e)  denote  the  right  hand  sides  of  formulas  (55),  (56), 
and  (57)  respectively.  Note  that  Ii{t)  and  /3(e)  can  also  be  expressed  as 


and 


or 


Thus 


fZTr 

h{t)=  I  fe{{x,ue) -\-t)d9 

h(e)=ri 

Vo  J\t 

f2Tr  r 


0  J\t\>e 

f°°  fe{{x,U9)  +  t) 


dtdO 


dtdO . 


=  h*f(x) 


and 


-1 


a,f(x)  =  g^{2/3(«)  -  2/2(e)} 


/„ 


1  /  1  v^|y/t|^  - 1 


x(y/^)  >  f{x  -  y)dy 


Ie?  \\yl€\  |2//e|2 

=  K,*f{x) 

which  are  the  desired  results. 

Formula  (56)  is  simply  a  consequence  of 

/•27r  r  2 

fe{{x,  ue))de  =  -^^/{x  -  y)dy 

and 

f  dt  2 
J\t\>e  e  ’ 

To  see  formula  (57)  use  identity  (55)  to  write 

2^  f  fe{{x,  Ue)  - f  f  f{x-y) 


Jo  j\t\>€ 


=  2  f  f 

J\t\>€  Jh 


l*l>^  4l>|t|  t^^J\y\‘^ 


dydt 
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=2/  (/ 

4i>«  (4 

=  2  I 

Jh 


dt 


|s/|>lt|>e  t'^^\y\^  —  ^2 

iyi>^  \yr 

which  is  the  desired  result. 

Finally,  to  see  (55)  first  recall,  from  Subsection  3.1,  that 


f{x  -  y)dy 


f27r  r27r  ^oo 

/  fei{x,  ue)  -  t)d9  =  /  /  f{x  -  tug  +  sv9)dsd9 . 

Jo  Jo  J—oo 


Now 


where 


-tug  +  svg  —  +  s^{a cos9  +  P sin 9,asm9  -  P cos 9) 


a  = 


v'PT' 


cos  (j)  and  P  = 


^/FT^ 


sin^ 


so  that 


-tUg  +  SVg  =  -Vt^  +  S^(cOs(9  —  (/>),  sin(^  -  ^))  —  +  S^Ug-4, . 


Thus  we  may  write 

c2'k  poo 


i'2'k  foo  roo  p2'7r  - 

/  /  f(x  —  tUg  +  SVg)dsd9  =  /  /  f{x  —  VW+  s'^ug-^)d9ds 

Jo  J—oo  J—oo  Jo 


poo  (  p2'K  _  ^ 

=  J  |y^  fix  -Vt‘^  +  s'^Ug)d9j 

~^lo  {lo 

27r  (roo  fdr  1 

f{x-rug)- 


p27r  (  poc 

[j\t\ 


y/r^  — 


ds, 

d9, 

d9, 


=  2  [ 
Jh 


fjx  -  y) 

'|y|>|t|  ^|i/|2  - 12 


dy 


where  the  second  equality  follows  from  translation  invariance  in  the  9  vari¬ 
able,  the  third  is  a  consequence  of  the  fact  that  the  integrand  is  even  in  s 
variable,  the  fourth  is  a  result  of  the  change  of  variable  r  =  y/t^  -I-  §2,  and  the 
last  follows  from  the  polar  change  of  variable  y  —  rug.  This  implies  (55)  and 
completes  the  proof  of  the  formulas  gef{x)  =  ke*f{x)  and  G^fix)  =  K(*fix). 

As  mentioned  earlier  the  various  assertions  of  the  Theorem  are  routine 
consequences  of  the  representation  G^f  =  K^*  f  and  the  properties  of  K. 
More  specifically: 
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•  Statements  ii,  iii,  and  iv  are  well  known  consequences  of  the  fact  that 
K  is  integrable  and  fj^2  K(x)dx  =  1. 

•  Statement  v  follows  from  the  additional  fact  that  /jjs  |a;|“|ii'(2:)|c?a;  is 

finite  if  0  <  O'  <  1  and  for  |a;|  >  2. 

•  Assertion  vi  follows  from  the  yet  additional  fact  that  A'  is  a  positive 
radial  function. 

•  Statement  vii  follows  from  the  general  theory  of  Holder  classes  of  func¬ 
tions  and  the  fact  that 

/  \Kix  +  y)-  2K(x)  +  K{x  -  y)\^dx  <  C\y\'^-^ 

JeP 

ii  1  <  q  <  2,  where  C  is  a  constant  independent  of  y.  For  example  see 
[2,  14,  24,  31]. 

•  Finally,  i  follows  from  the  facts  that  Jjf  K{x)dx  =  1  and  that  |A'(a:)|  is 
dominated  by  a  integrable  radially  decreasing  function,  see  [31].  Note 
that  k{x)  does  not  enjoy  the  second  property. 

3.4 

There  are  many  ways  of  arriving  at  (16).  Perhaps  the  most  direct  and  sim¬ 
plest  method  is  to  use  (i)  the  fact  that  the  (univariate)  Fourier  transform 
of  fe  is  equal  to  the  (bivariate)  Fourier  transform  of  /  restricted  to  the  line 
TUe,  -oo  <T  <00,01  symbolically  /^(r)  =  firue),  (ii)  the  Fourier  inversion 
formula,  and  (iii)  a  polar  change  of  variables.  For  example  see  [22,  29].  It  is 
also  a  natural  consequence  of  various  summability  formulas  for  the  inversion 
of  the  Radon  transform,  see  [20,  page  196].  Of  course  it  can  also  be  derived 
from  (1)  by  using  the  fact  that 

/e(t)  = 

and  an  interchange  of  order  of  differentiation  and  integration. 

Reasonable  smoothness  conditions  on  /  needed  to  make  formula  (16)  valid 
are  not  clear.  Certainly  /  6  for  some  positive  a  is  suflicient^^  but  is 
far  too  restrictive  for  typical  applications.  For  this  reason  we  take  (19)  to  be 
the  definition  of  Kf{x)  and  view  (18)  as  an  expression  valid  for  sufficiently 
smooth  /. 

11(71+“  is  the  class  of  those  differentiable  functions  whose  derivatives  are  Holder  con¬ 
tinuous  of  order  a. 
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3.5 


Relationship  (22)  is  valid  whenever  ^  is  a  reasonably  smooth  radial  function, 
it  need  not  necessarily  be  integrable.  Indeed,  in  [20]  this  relationship  is 
systematically  exploited  with  polynomial  ^’s. 

On  the  other  hand,  (24)  need  not  even  make  sense  if  ^  is  not  integrable. 
Nevertheless  if  ^  and  sufficiently  many  of  its  derivatives  are  integrable  then 
not  only  is  xp  well  defined  by  (24)  but  it  is  also  integrable  and  has  mean  value 
zero.  In  this  case,  even  if  ^  fails  to  be  an  analyzing  wavelet  because  ^(0)  ^  0, 
the  corresponding  xp  is  always  a  univariate  analyzing  wavelet.  Thus  if  is 
a  well  behaved  approximation  of  the  identity  (23)  implies  that  '^a  *  /(^)  can 
be  computed  as  an  average  of  wavelet  transforms  of  /$,  0  <  6  <  tt,  and 
consequently  (23)  may  be  regarded  as  an  approximate  inversion  formula  in 
terms  of  wavelet  transforms. 

3.6 

The  way  (26)  is  computed  in  terms  of  the  Radon  transform  data  of  /  is 
dictated  by  formula  (8)  and 

($,*/),(t)  =  ($,)e*/fi(t). 

This  results  in 

A(«,  •  /)(^)  =  -^  •  /((*.  ■ 

For  more  details  see  [7,  13,  30].  Thus 


Mt)  = 


which  is  another  variant  of  (29). 


3.7 

To  see  the  formula  for  xp  observe  that  for  any  sufficiently  well  behaved  uni¬ 
variate  function  g 

=  -/  a(s  +  t)-2g(s)  +  gis-t)^^ 

Ja/ 


TT 


gjs  - 1) 

i/2<t<o  a?{t/aY 

-ypa*  9{s) 


477  J a/2<t<.a 

gjs)  a 

TT 

p(s) 


-/ 

27r  Ja/ 


dt 


30 


where 


i)(t)  = 


xi2t)  -  x{t) 

27ri2 


■00 (t)  =  a~^'ip{t/a),  and  x(2i)  —  x(^)  is  the  indicator  function  of  {f  :  1/2  < 
1^1  <  !}• 


3.8 


There  are  several  points  of  view  concerning  convolution-backprojection  type 
reconstruction  algorithms  in  computed  tomography,  [4,  6,  10,  13,  20,  22,  23, 
27,  28,  29,  30].  We  will  not  outline  them  all  here  and  only  mention  that  the 
connection  between  the  various  discrete  implementations  and  the  inversion 
methods  on  which  they  are  based  have,  to  my  knowledge,  not  been  defini¬ 
tively  developed.  These  connections  are  usually  based  on  well  established, 
but  nevertheless  heuristic,  principles  rather  than  precise  error  estimates.  The 
relatively  recent  papers  [6,  23]  nicely  illustrate  the  state  of  the  art. 

For  example,  the  artifacts  introduced  by  discretization  in  the  9  variable  of 
the  commonly  used  inversion  methods  are  well  known.  However,  they  have 
not  been  quantitatively  analyzed.  On  the  other  hand,  the  reconstruction 
algorithms  introduced  in  [15, 17]  do  not  give  rise  to  such  artifacts,  indeed  this 
is  one  of  the  central  points  of  that  work.  Unfortunately,  unlike  some  of  the 
other  methods  described  here,  these  algorithms  have  no  localized  analogues. 

Since  [16]  may  not  be  easily  accessible,  we  now  briefly  outline  the  rationale 
for  the  paradigm  used  to  approximate 

/(X> 

w{s  —  t)g{t)dt 

-OO 

where  w  is  a  known  univariate  integrable  function  or  measure,  whose  deriva¬ 
tive  is  very  large  or  even  non-existant  in  the  classical  sense,  and  is  a  uni¬ 
variate  function  supported  in  the  interval  —  1  <  t  <  1  whose  samples  g{n/N), 
n  =  0,  ±1,  ±2, . . . ,  ±N,  are  known. 

First  assume  that  g  is  reasonably  well  behaved  so  that  approximation 

12  g{n/N)X{t-n/N) 

n=-N 

enjoys  the  correct  order  of  approximation.  For  example,  if  g  is  Holder  con¬ 
tinuous  of  order  a  then 

(58)  \9{t)-m\<j;^- 

The  primary  rationale  of  this  method  is  that  the  approximant  be  easy  to 
compute  and  that  the  approximation  of  w  *  ^  be  of  the  same  order  as  \g{t)  — 
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g{t)\,  e.  g.  (58).  Now,  consider  the  method  outlined  in  Subsection  2.8,  that 
is,  replace  w  with 

CX) 

w{t)=  Y, 

n=~cxD 

where  6{t)  is  the  univariate  unit  Dirac  delta  “function”  at  the  origin, 

/OO 

w{t)(l)n(t)dt , 

-OO 

and  <l)n{t)  is  sufficiently  smooth,  is  centered  and  concentrated  around  n/N, 
and  has  total  integral  one.  It  is  not  difficult  to  see  that  (i)  w  *  g  is  of  the 
same  form  as  g,  that  is,  it  is  a  linear  combination  of  translates  of  A,  (ii)  if 
the  (j)nS  are  chosen  appropriately  then  w  *  g  is  easy  to  compute,  and  (iii) 
\w  *  g{t)  —  w  *  g{t)\  is  of  the  same  order  as  \g{t)  —  g(t)\,  e.  g.  if  (58)  is  the 
case  then 

\w  *  g(t)  -  w  *  g{t)\  < 
where  L  is  the  norm  or  total  variation  of  w. 
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Figure  1:  The  reconstruction  Af{x)  via  formula  (43). 
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Figure  6:  Pull  data,  conventional  reconstruction,  and  local  wavelet  transform 


Figure  7:  3/4  of  full  data,  conventional  reconstruction,  and  local  wavelet 
transform.  Scaled  to  occupy  the  same  area  as  Figure  6. 


Figure  8:  1/2  of  full  data,  conventional  reconstruction,  and  local  wavelet 
transform.  Scaled  to  occupy  the  same  area  as  Figure  6. 


37 


Figure  9:  1/4  of  full  data,  conventional  reconstruction,  and  local  wavelet 
transform.  Scaled  to  be  proportional  to  Figure  8. 


Figure  10:  1/8  of  full  data,  conventional  reconstruction,  and  local  wavelet 
transform.  Scaled  to  be  proportional  to  Figures  8  and  9. 
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Figure  11:  Cross  sections  corresponding  to  the  reconstruction  in  Figure  6 


38 


